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IN MEMORY OF RYZHKIN 





On January 27, 2021, at the 84th year of his life, a remarkable scientist, engineer and 
organizer, Doctor of Engineering Sciences, Professor Anatoly Andreevich Ryzhkin died a 
sudden death. 

Anatoly A. Ryzhkin was Honored Worker of Science and Technology of the Russian 
Federation, Doctor of Engineering Sciences, Professor. 

Anatoly Andreevich was born on January 21, 1938. He graduated from Rostov Institute 
of Agricultural Engineering (RISKHM) in 1960. All of his more than 60-year working 
career has been associated with his native university. 

In the transition time for RISKHM and the entire Russian education, he headed the 
University — in February 1988, he was elected Rector and held this position until 2007. 
Anatoly A. Ryzhkin made a great contribution to the transformation of a branch institute of 





agricultural engineering into a technical university: a modern base for scientific and 
educational process was created, the number of specialties was increased, postgraduate programs were expanded, and 
doctoral programs were opened, social infrastructure facilities were developed. As a result, RISKHM has become a 
regional center of education, science, and culture. In 1992, on his initiative and under his personal leadership, by the 
decision of the Government of the Russian Federation, RISKHM was transformed into Don State Technical University 
(DSTU). 

Anatoly A. Ryzhkin is a well-known scientist in the field of friction and wear under the conditions of cutting 
materials. He founded and headed the scientific direction “Physicochemical and thermodynamic fundamentals of wear 
resistance control of tool cutting materials”. A. A. Ryzhkin's applied research works were aimed at creating 
technologies for strengthening tool materials. The results of R&D are implemented at various enterprises of the Rostov 
region and are widely used in the training of students majoring in “Design and technological support of machine- 
building industries”. 

A. A. Ryzhkin is the author of more than 450 scientific and method papers. More than 50 of them were 
published abroad in English, German, Polish, Spanish, and other languages. Over the years, Professor A. A. Ryzhkin 
had led the postgraduate and doctoral studies, was Chairman of the Dissertation Defence Board at Don State Technical 
University. Anatoly Andreevich personally prepared 16 candidates and 5 doctors of engineering sciences. 

For more than 15 years, Anatoly A. Ryzhkin headed the Council of Rectors of Universities of the Rostov 
region, which, under his leadership, became an effective public and state body for governing the education sector of the 
Rostov Region. 

In recent years, Anatoly Andreevich was Head of the Metal-Cutting Machines and Tools Department, as well 
as the Vice-president of the Association for Engineering Education of Russia, a member of the Problem Council for 
Mechanical Engineering under the Ministry of Education and Science of the Russian Federation. 

For long-term and conscientious work, a great contribution to solving problems of higher education, large- 
scale scientific and public activities, Anatoly A. Ryzhkin was awarded various state awards and honorary titles: Veteran 
of Labor, Honorary Professor of DSTU, Honored Worker of Science and Technology of the Russian Federation, 
Honorary Worker of Higher School of the Russian Federation. He is laureate of awards of President of the Russian 
Federation and the Government of the Russian Federation in the field of education; holder of the Order of Merit for the 
Fatherland, IV degree, the Order of Friendship, the Badge of Honor, the Medal for Valiant Labor, the Order of Merit to 
the Rostov Region. 
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Research of physicomechanical and design characteristics of vibrated, 
centrifuged and vibro-centrifuged concretes 
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Introduction. Currently, the obtaining of lightweight concrete and reinforced concrete products and structures with the 
improved structure and characteristics is a challenge. This can be achieved through centrifugation or in a more 
advanced way — vibro-centrifugation. At the same time, the influence of centrifugal and centripetal forces of inertia in 
these types of technologies causes differences in the cross-section properties of concrete products and structures. To 
reflect this in the calculations, it 1s required to experimentally and analytically investigate the qualitative and 
quantitative patterns of such differences in the characteristics of concretes obtained through different technologies. 
Materials and Methods. The study used the cross-section averaged characteristics of concrete — “integral 
characteristics of concrete”. The applicable raw materials included portland cement 500, crushed stone fraction 5-20, 
medium sand. Nine control samples of annular cross-section obtained through vibrating, centrifuging, and vibro- 
centrifugation were manufactured and tested. The essence of the technique was that each manufactured experimental 
control sample was used in several types of tests in-parallel. From the total annular section of each sample, three 
conditional quadrants were distinguished, from which standard samples of small size were cut out. Subsequently, they 
were tested for axial compression, tension, and flexural tension. The following test equipment was used: electronically 
controlled mechanical press IPS-10 — for compression testing of prisms, and the breaking machine R-10 — for testing 
samples for axial tension. Strain sensors and dial indicators were used to measure concrete deformations. Oscilloscopes 
were also used to obtain the deformative and strength properties of concrete, including full deformation diagrams with 
descending branches. 

Results. We have analyzed the calculation results of the integral design characteristics of the concretes obtained through 
vibration, centrifugation and vibro-centrifugation. It is established that due to the influence of centrifugal and centripetal 
forces of inertia under centrifugation and vibration centrifugation, the characteristics of concrete in cross-section 
become different. In some cases, these differences can be very significant. We have developed and tested the following: 
a new method for evaluating the dependence of the integral (cross—section averaged) design characteristics of concrete 
(density, cubic and prismatic axial compressive strength); ultimate deformations under axial compression; axial tensile 
and flexural tensile strength; ultimate deformations under axial tension; elasticity modulus; diagram of “stress 6,— 


7 


strain €,” under compression; diagram of “stress 6,,—strain &, under tension on the manufacturing technology 
(vibrating, centrifuging, vibration centrifugation). 

Discussion and Conclusions. Based on the results of the research, conclusions are formulated on the positive effect of 
the proposed technology of joint vibrating and centrifuging. It consists in improving the integral design characteristics 


and structure of concrete from vibrating to centrifuging and from centrifuging to vibro-centrifuging. 


Keywords: vibrating, centrifugation, vibro-centrifugation, column calculation, variatropic structure, integral 


characteristics of concrete, ultimate deformations, compressive strength, elasticity modulus. 


For citation: L. R. Mailyan, S. A. Stel'makh, E. M. Shcherban', et al. Researches of physicomechanical and design 


characteristics of vibrated, centrifuged and vibro-centrifuged concretes. Advanced Engineering Research, 2021, 
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Introduction. It is known that through vibro-centrifugation'’”” [1-11], it is possible to obtain concretes with 
improved structure and characteristics compared to those obtained through centrifugation and vibration. 

However, in a few works on concrete and reinforced concrete annular structures, the influence of the type of 
technology on the average (general) cross-section characteristics of concrete was considered [12—15]. At the same time, 
it is obvious that due to the influence of centrifugal and centripetal forces of inertia under centrifugation and vibration 
centrifugation, the cross-sectional characteristics of concrete differ [16]. 

In this paper, we experimentally and analytically investigate the qualitative and quantitative patterns of such 
differences in the characteristics of concretes obtained through various technologies. Obviously, to account for these 
differences, it is required to introduce some averaged cross-sectional characteristics of the elements. For this purpose, 
we introduce the term “integral characteristics of concrete”’. 

Materials and Methods. In total, nine control annular samples made through vibration, centrifugation and vibro- 
centrifugation were produced and tested. Dimensions of these samples were as follows: outer diameter D — 450 mm, 
inner diameter d — 150 mm; total height H — 1200 mm. 

The equipment and test methods used are described in [8-11]. 

Crushed stone fraction 5—20 was used as a filler, which brings the properties of the resulting concrete closer to 
the properties of fine-grained concrete. 

In the experiments, the type of manufacturing technology varied, which was recorded in the designations of the 
samples: vibration — V, centrifugation — C, vibro-centrifugation —VC. 

The problem of estimating the dependence of the integral (averaged over the cross-section) structural 
characteristics of concrete (density, axial compressive strength (cube and prism); ultimate deformations under axial 
compression; axial tension strength and flexural tension; ultimate deformations under axial tension; elasticity modulus; 
diagram of “stress o,— strain ¢,’ under compression; diagram of “stress G,—strain ¢€,,’ under tension) on the 
manufacturing technology (vibration, centrifugation, vibro-centrifugation) was considered. 

Research Results. The test procedure differed in that each manufactured experimental control sample was 
used in several types of tests in-parallel. Control samples in a single copy were selected and tested on the 7th, 28th and 
180th days. 

From the total annular section of each sample, 3 conditional quadrants A, B and C were allocated, from which 
small-size samples were cut out. Subsequently, they were tested for axial compression, tension, and flexural tension 
(Fig. 1, 2). 

Four cube samples with an edge of 15 cm were cut out of the quadrant A for compression and tension tests 
(levels 1-4), one prism (15<x15x60 cm) — for flexural tension tests (level 5). For the axial compression tests, two 
prisms (15<x15x60 cm) were cut out of the quadrant B (levels 1—2). Next, three prisms (15<15x60 cm) were cut out of 
the quadrant C for axial tension tests (levels 1—2). 

After testing the cubes for axial compression, we obtained values R,.,,,, after testing the prisms for axial 
compression — values Rp, Epp, Ry, Exir, E;=Ep, and the deformation diagram “6,—€,”; after testing the prisms for axial 
tension — values R,, and the deformation diagram “6,,;-€;;’; and after testing the prisms for flexural tension — values 


Roz. 


' Aksomitas GA. Strength of short centrifuged annular section columns with longitudinal reinforcement of class At-V under short-term compression: 
Cand.Sci. (Eng.) diss. Vilnius: VISI; 1984. 261 p. (In Russ.) 

* Petrov VP. Technology and properties of centrifuged concrete with combined aggregate for overhead contact line supports: Cand.Sci. (Eng.) diss. 
Rostov-on-Don: RISI; 1983. 175 p. (in Russ.) 

> Rajan Suwal. Properties of centrifuged concrete and improvement of the design of centrifuged reinforced concrete poles for power transmission 
lines: Cand.Sci. (Eng.) diss. Rostov-on-Don: RGAS; 1997. 267 p. (In Russ.) 
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Fig. 1. Scheme of manufacturing small-size concrete samples from quadrants along the height of experimental control full-scale 
samples of annular cross-section for calculating integral characteristics: a) quadrant A; b) quadrant B; c) quadrant C 
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Fig. 2. Experimental small-size concrete samples from quadrants A, B, C of experimental control samples of full-scale annular cross- 
section for the analysis of integral structural characteristics a) quadrant A; 6) quadrant B; c) quadrant C 


The test procedure according to GOST 10 180 was applied. The following test equipment was used: IPS-10 — 
for testing prisms for compression, and P-10 — for testing samples for axial tension. 

To measure concrete deformations, strain sensors with a measurement base of 50 mm and dial indicators with 
graduation 0.001 mm were used. 

To obtain the deformative and strength properties of concrete, including complete deformation diagrams with 
descending branches, experiments were carried out using a constant deformation rate. 

For this purpose, in addition to strain gauges, oscilloscopes were also used. 

For tests with the same load feed rate, a loading step of 0.1R was selected, and the deformations of the prisms 


increased with a step of 0.1€, (Fig. 3). 
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Fig. 3. Test modes of experimental small-size samples: 
a) — step-by-step, with constant loading rate (Ac=const); 
b) — step-by-step, with constant deformation rate (Ae=const) 


The final test mode consisted in increasing the load to the maximum and its subsequent decrease in the process 
of increasing deformation. Thus, during the tests, the descending branch of the concrete diagrams “6 — ¢” was fixed, 
which has a fairly clear outline up to about the value of o = 0.8 R, both under compression and tension. Subsequently, 
the dependence acquired a very unstable character 

After analyzing the results obtained, we can draw conclusions about the influence of the sample manufacturing 
technology. The results of experimental studies of changes in the integral characteristics of experimental concrete 
samples, depending on the manufacturing technology, are shown in Fig. 4—10. 
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Fig. 4. Influence of manufacturing technique on increase in cube strength of concrete under compression 
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Fig. 5. Influence of manufacturing technology on increase in prism strength of concrete under axial compression 
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Fig. 6. Dependence of increase in bending tensile strength of concrete on manufacturing technique 
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Fig. 7. Dependence of increase in axial tensile strength of concrete on manufacturing technique 
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Fig. 8. Dependence of reduction of ultimate deformations under axial compression of concrete on manufacturing technique 
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Fig. 9. Dependence of reduction of ultimate deformationsunder axial tension of concrete on manufacturing technique 
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Fig. 10. Influence of manufacturing technique on increase in the concrete modulus of elasticity 


Discussion and Conclusions. The influence of the sample manufacturing technology on the density at all ages 
was minimal (about 2 %), therefore, when calculating, the indicator “density” can be neglected. 

Vibro-centrifugated concretes showed higher values in terms of “compressive strength” and “tensile strength” 
than vibrated and centrifuged concretes, namely: prism and cube compressive strength — up to 22.0 %, axial and 
flexural tensile strength — up to 27.0%. 

Due to the ordering of the ongoing processes of hydration of cement stone, an increase in the compressive and 
tensile strengths of concrete with simultaneous vibration and centrifugation, 1n comparison with centrifuged and 
vibrated samples, with age occurs in the studied range of concrete age (7-180 days) and is practically within the same 
limits in all ages. 

For vibro-centrifuged samples, there is a slight decrease (up to 6 %) in the ultimate deformations under axial 
compression and tension. This corresponds to the maximum strength of the concretes under study. 

Vibro-centrifugated concretes, in comparison with vibrated and centrifuged concretes, showed the lowest 


ultimate deformations at any age. 
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The values of “compressive modulus E,” and “tensile modulus E,,” at all ages of concrete were up to 8.0 % 
higher in vibro-centrifuged concretes than in vibro-centrifuged and centrifuged concretes. 

An increase in strength, with a parallel decrease in the ultimate deformations, was the reason for an increase in 
E, and E,,1n concretes with simultaneous vibration and centrifugation, rather than in concretes with one of the types of 
compaction. This fact affected the “stress-strain” diagram — the maximum shifted up and to the left. 

Differences in the stress-strain diagrams under compression and tension, characteristic of concretes with 
simultaneous centrifugation and vibration, at all ages are as follows: 

— increase in strength and decrease in ultimate deformations (maximum shifts up and to the left); 

— increase in the initial modulus of elasticity (increase at the ascending angle origin). 

At all ages of concretes, the following tendency was characteristic: an increase in the lifting capacity of the 
ascending branch of the diagrams, a decrease in the descending branch in centrifuged and vibro-centrifuged concretes 
compared to vibrated ones. 

According to this study, it is appropriate to draw the following conclusions. 

1. Studies of the integral design characteristics under compression and tension of the considered types of 
concrete of various manufacturing techniques at the ages of 7, 28 and 180 days showed: 

— concrete performance is improved from vibrating to centrifuging and from centrifuging to vibro- 
centrifuging; 

— increase at all ages of compressive and tensile strengths (up to 23 %); 

— reduction of all ultimate deformations (up to 8 %); 

— increase in elastic modulus under various types of loading (up to 10 %). 

2. For all integral diagrams of the “stress-strain” deformation of concrete with simultaneous vibration and 
centrifugation, the following is characteristic: 

— moving the maximum up and to the left; 

— increasing at the ascending angle origin; 

— increasing the chart lift in the ascending branch. 

3. Through numerous experimental studies, it has been established that simultaneous vibration and 
centrifugation contributes to the production of concretes with improved structure and characteristics, rather than 
concretes obtained by only one type of impact — centrifugation or on. 
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Three-dimensional integral dry friction model for the motion 
of a rectangular body 





M. S. Salimovy, I. V. Merkuriev 
National Research University "MPEI" (Moscow, Russian Federation) 


Introduction. A three-dimensional dry friction model in the interaction of a rectangular body and a horizontal rough 
surface is considered. It is assumed that there is no separation of the body from the horizontal surface. The body motion 
occurs under the conditions of combined dynamics when, in addition to the longitudinal movement, the body partici- 
pates in twisting. 

Materials and Methods. Linear fractional Pade approximations are proposed, which replaced the cumbersome analytical 
expressions that most accurately describe the motion of bodies on rough surfaces. New mathematical models describing 
sliding and twisting of bodies with a rectangular base are proposed. 

Results. Analytical expressions of the principal vector and moment of friction for rectangular contact areas are devel- 
oped and scientifically established. A friction model that takes into account the relationship between sliding and twist- 
ing speeds, which provides finding solutions for Pade dependences, is developed. After numerical solution to the equa- 
tions of motion, the dependences of the sliding speed and angular velocity on time were obtained and constructed. 
Graphs of the dependences of the friction forces and their moment on two parameters (angular velocity and slip veloci- 
ty) were constructed, which enabled to compare the integral and normalized models of friction. The comparison results 
showed good agreement of the integral model and the model based on Pade approximations. 

Discussion and Conclusions. The results obtained provide considering the dynamic coupling of components, which 
determines the force interaction of a rectangular body and a horizontal surface. These results can be used in mobile ro- 
botics. The analyzed motion of the body occurs through the motion control of a material point inside the body. Such 
mobile robots can be used when solving a wide class of problems: when creating autonomous robots for the exploration 
of outer space and planets; in the diagnosis and treatment in case of passing through complex structures of veins and 
arteries; in research under water, in places of large differential temperature; in underground operations. 
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Introduction. The study of the movement of a rectangular body is a challenge in the mobile robotics [1]. This 
movement is due to the control of the material point inside the body. Such mobile robots can be used to solve a wide 
range of tasks. For example, when creating autonomous robots for the exploration of outer space and planets; for 
medical purposes, in diagnosis and treatment, for example, in case of passing through complex structures of veins and 
arteries; as well as for underground work and research under difficult conditions, for example, under water and in places 
of large differential temperature [1, 2]. 

Thus, more and more challenges are being set for robotics, which require theoretical research, including 
studying models of friction between the body and the surface under the conditions of combined dynamics [3, 4]. Since 
the movement of the mobile robot occurs in different directions, it is required to consider the longitudinal movement 
and rotation. Thus, in the structure of the friction model, it is required to provide the relationship between the sliding 


and twisting speeds [5]. An important development in the description of this relationship was made in [6]. Its author 
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managed to solve the equations for the principal moment and the vector of friction forces where a rectangle was 
considered as a contact area. Such analytical expressions enable to most accurately describe the motion of bodies on 
rough surfaces, but they are cumbersome and complex since they contain integral expressions. Hence, the authors of [7] 
constructed linear fractional Pade approximations, which made it possible to find solutions for the resulting 
dependences. 

Pade approximation can be used to explain the effects of combined dry friction for linear and angular 
velocities. On the basis of Pade approximations, it became possible to create new models of friction [8, 9], which later 
began to be classified for better interpretation [10]. The classification occurs depending on the number of parameters. 
Thus, in [11], the authors introduced the notions of dimension and order of the dry friction model depending on the 
order of the used Pade approximations. 

The model of sliding and twisting friction, which is proposed in the paper [12], provides considering the 
dynamic connection of the components that determine the force interaction of a rectangular body and a horizontal 


surface [13]. 


Problem Statement. We consider a solid body of mass m, , which is a rectangular body with uniform faces of 
length a, width 5 and height 2h. A fixed coordinate system Oxyz , associated with the body (Fig. 1) is introduced. 
Point O is located on the horizontal plane. The system. O,x,y,z, starts at a point O,, that corresponds to the geometric 
center of the body. The axis Oz, is parallel to the axis Oz. The axis O,x, 1s parallel to the long edge of the body. We 


introduce the unit vectors e.,e, of the axes O,y, and O,x,, respectively. 
were 








Fig. 1. Coordinate systems 


Consider the continuous motion of the body on the surface (Fig. 2), which consists of translational movement 
and rotation about the axis O,z,. Three coordinates determine the position of the body. The coordinates x,, y, and A 


set the origin of the coordinate system O,x,y,z, in the coordinates Oxyz . The rotation of the body relative to its initial 
position on the axis O,x, is specified by the angle @. This paper considers the case when the center of mass of the body 


G and the center of mass of the system O, coincide (Fig. 2) [14]. 


2) 


xy 





mg 


Fig. 2. Movement of the system body 
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Materials and Methods. The contact area is a rectangle with sides a and b, in which the normal voltage 
depends on the distance from the point P to the faces of the rectangle (Fig. 3). 





Fig. 3. Velocities of the points O, and P 


Consider an infinitesimal area dS at an arbitrary point M on the contact surface. We introduce the angle 
between the relative sliding velocity and the axis O,x,. Let us draw the radius vector 7,,, from the point P to the point 


M. The velocity vector of the point M is denoted uy, and to find it, we use the Euler formula describing the velocity 
distribution in a perfectly rigid body: 
Uy = Up + W-Typ. 
The sliding speed at the point M is decomposed into two components along the axes 0,x, and O,y,: 
Umx = Vx — YO; 
Umy = Vy + XW. 
Using Coulomb's law, we find a small increment of the friction force directed against the relative velocity at 
the point M [15]: 
dF = —fo(x, yids, 


where f — coefficient of friction; o(x, y) — contact stress distribution function depending on the x and y coordinates; 


dS = dxdy — small area increment [15]. 


We rewrite the differential of the friction force and the moment of this force in projections on the axes under 





consideration: 
dE, = —fo(x, y)—™ dxdy; 
luv 
dF, = —fo(x,y) -“Y dxdy: 
lum 
i jk 
dM,=|*xX yy O|=xdF, — ydF,. 
dk, dF, 0 


As a special case, we consider a uniform distribution of stress in the absence of internal masses in the body, 
then these stresses will be equal to: o = a but then we will continue the record in general form: o(x, y). 


Having integrated the expressions for the friction forces, we obtain: 








/2 pb/2 = 
F, = —f ‘en | 5/3 o(x,y) ar dxdy; (1) 
a/2 -b/2 Vytxw 
Fy a = ine Spyz o(x,y) i dxdy. (2) 


Relative slip module |v,,| is calculated from the formula: 


luyl = v3, + Viry = v2 gO (ear y-) 2w(v,x — Uxy). (3) 
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Imagine the relative positions of the vectors of variable sliding speed uv and the components of the friction 
force: F — the component opposite to the sliding speed v; F,; — the component perpendicular to the instantaneous slip 


velocity. At the same time, imagine the coordination of this system with respect to the axes O,x, and O,y, (Fig. 4). 





Fig. 4. Components of the friction force and velocity 


We will make a transition from the projections of the sliding speed: 


Vy = UCOS®, 
is = vsing, (4) 


to the speed module and the sliding angle: 


F, = F,cos@ + F,sing, 


F, = F,(-sing) + Fsing’ (5) 


We will integrate the moment of the friction force on the contact area: 
a/2 (gee u(xsing—ycos@)+w(x2+y7) 


M, = —f Jays —p/2 a(x, y) xdy. 


v2+02(x2+y2)+2wvu(xsing—ycos@) 


Let us substitute the expressions (1)—(3) into the system (5), and also rewrite the expression for the moment of 
force. As a result, we obtain a three-dimensional model of friction sliding and twisting: 
a/2 


b/2 v(cos?@+sin? ~)—w(ycosptxsing) 


Py cae’ d lea ye o(x, y) v24+w2(x2+y2)+2wv(xsin~—ycos@) (6) 
= a/2 b/2 w(ysin@+xcosq@) 

P= poe? Son72 o(x, y) v24+w2(x2+y2)4+2wv(xsinp—ycos@) 7) 
_ _ ¢ fa/2 pb/2 u(xsing—ycos@)+w(x2+y7) 

M, = Pye Son72 o(x, y) v24+w2(x2+y2)+2wv(xsing—ycos@) aay co 


In order not to solve cumbersome integrals, we use the replacement of the corresponding Pade expansions [16, 
17]. Thus, based on the Pade theory [18], these expressions can be formulated as the ratio of two functions of several 
variables in the entire domain of definition, provided that the functions must have the same order [7]. To define these 


functions, it is required to determine the behavior of integral expressions (6)—(8) under the following conditions: 











OF | _ _ f a/2 pb/2 x2 +y?(x*+y7)+(ycosptxsing) (xsinp—ycos@) =: oy. 
dv \v=0 7 w ?—-a/2 ee a(x, y) (x2+4y2)2 dxdy ~ Ww Io; 
aM, _ _ f pa/2 b/2 _ 2y?xsing—2y* cos@) a ae 
dv lu=0 wt -a/2 pe o(x, y) (x2 4+y2)2 dxdy 133 
a/2 ¢b/2 x2 —y? 

Mo e900 a =f Joss eee o(x,y) " Were iad = = fle. 

a/2 b/2 ysin~+xcos@ 

Figoo = —F Saj2 Sp). TY) Tlgye  OXAY = flo; 
oh set ser 2 rr ae 
OV jo=0 IW|y=0 IW|w=0 Nese OV Jw=0 IW |y=0 

_ _ OF, _ OM, _ OM, _ _ 

—_ Pies ~ da l=0 ~ Oy io=6 Ba lv=0 _ M7 1y-00 = 0. 
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. OF aM 
Values of the expressions —= = 


are not involved in finding the subsequent Pade approximants, 
dV |v=0 dW |w=0 





therefore, their writing is omitted due to their cumbersomeness. The identical equality to zero is realized under the 
condition that the voltage o is symmetric about the center of the rectangular contact spot, 1.e., the point P. 

An accurate three-dimensional integral model [13] (6)-(8) provides a logical description of dry friction 
phenomena, but for solving problems of dynamics, such a model is difficult to accept due to the need to calculate 
impressive integrals [10]. To avoid this procedure, we use [6] to replace the exact integral system with the 
corresponding expressions using Pade approximations in the whole range of variables. Linear fractional Pade 


expansions give a three-dimensional model of first-order sliding and twisting friction [19]: 











v+b1w 
Fi = Fo ao (9) 
_ W+b2v_ 
M, = Mo (10) 
_ W+b3u 
Td aera (11) 


To determine the Pade coefficients, it is required to study the properties of this model at the boundary points 
by analogy with integral expressions. To do this, we differentiate the parameters F), F,, M, and thus satisfy the 


corresponding integral expressions: 


(Fai —— 


fo 
vt T, WO 
M, = —flg 
I OO 
Fy = -fle—— 
{ w+ Fv 
9 
The system of equations of motion has the form: 
aWo, 
J dt Z 
dV x 
(mp +m) _ = F, + (mp + My, )VyWo,; (12) 
dvy 
(Mm, + m,) “ae _ ky —- (mM, at M)Vx.Wo,- 
We express the time derivatives of the sliding speed and the angular speed using the formulas (3)-(5): 
du 1 ; dv, 5 dv, 
— = ———— | 20, —— + 20, — |; 
dt 2./v2 +02 ~ dt Y dt 
d 1 v sin v cos 
a (Bain ens 
dt m)j+m, vu? v2 


We rewrite these equations using the formulas (3)-(5) for Wo, = w/a and add the first equation from the 
system (12): 
J ad = M,a; 
dt — 
(Mp + my) < = Fy; (13) 
(m,) + m,)Uuq@ = F,. 


Research Results. Next, we calculate the integral expressions of the parameters Jp, 1,, Iz, Ig, Ig c using the 


Wolfram Mathematica software package for the following values: 


f=1,a= 05m; b= 0.2m; my = 1kg, o= a a ee 
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and substitute in the system of equations (13). Based on numerical expressions, we build graphs of integral and 
normalized functions depending on the parameter = ¥/,, . Fig. 5 shows graphs of the functions of the integral friction 


els (11)-(13), as well as models based on Pade approximations (14)-(16). 


Fy Mc 
0,6 
0.020 
0.5 Fh 


0.015 | : 
0.010 


0.005 

















a) k b )) k 
Fig. 5. Graphs of integral (solid lines) and normalized (dotted lines) functions of the tangent (a), normal (b) components of friction 
force and friction moment (c) 


Based on the graphs of the functions (Fig. 5), we can talk about good matching of the considered models. 


Next, we obtain graphs of the dependences of the characteristic parameters on time (Fig. 6). 





v(2), m/s w(t), m/s 

0.20 ° 
4) 
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0.05 1 | 
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Fig. 6. Dependences of sliding velocity v (a) and angular velocity w (b) on time ¢ 


Discussion and Conclusions. The movement of the mobile robot, starting from the contact of its body and the 
reference plane, under the conditions of combined dynamics, when there is sliding and twisting, is described. Analytical 
integral expressions are obtained for the tangent and normal components of the friction force [19] and the moment of 
friction applied to a rectangular contact area. The corresponding Pade approximations are determined for the obtained 
expressions. The integral and normalized models are compared through plotting the dependences of the friction forces 
and the moment of friction on the angular velocity and the slip velocity. The comparison results showed good matching 
of the integral model and the model based on the Pade decompositions. The graphs correspond to the logical behavior 
when a rectangular body moves, since the sliding speed and angular velocity increase according to the specified 
parameters. Consequently, the combined friction model implemented using Pade approximations can be applied to 
solve problems related to mobile robots with a rectangular base. 
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Thick-walled spherical shell problem 
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Introduction. Cylindrical and spherical shells are extensively used in engineering. They face internal and/or external 
pressure and heat. Stresses and strains distribution in elastoplastic shells has been studied by many scientists. Numerous 
works involve the use of the von Mises yield conditions, maximum shear stress, maximum reduced stress. These condi- 
tions do not include the dependence on the first invariant of the stress tensor and the sign of the third invariant of the 
stress deviator. In some cases, it is possible to obtain numerical-analytical solutions for stresses, displacements and de- 
formations for bodies with spherical and cylindrical symmetry under axisymmetric thermal and force action. 

Materials and Methods. The problem on the state of a thick-walled elastoplastic shell is solved within the framework of 
the theory of small deformations. A plasticity condition is proposed, which takes into account the dependence of the 
stress tensor on three independent invariants, and also considers the sign of the third invariant of the stress deviator and 
translational hardening of the material. A disconnected thermoelastoplastic problem is being solved. To estimate the 
stresses in the region of the elastic state of a spherical shell, an equivalent stress is introduced, which is similar to the 
selected plasticity function. The construction of the stress vector hodograph is used as a method for verification of the 
Stress state. 

Results. The problem has an analytical solution for linear plasticity functions. A solution is obtained when the strength- 
ening of the material is taken into account. Analytical and graphical relationships between the parameters of external 
action for the elastic or elastoplastic states of the sphere are determined. For a combined load, variants are possible 
when the plastic region is generated at the inner and outer boundaries of the sphere or between these boundaries. 
Discussion and Conclusions. The calculation results have shown that taking into account the plastic compressibility and 
the dependence of the plastic limit on temperature can have a significant impact on the stress and strain state of a hollow 
sphere. In this case, taking into account the first invariant of the stress tensor under the plasticity condition leads to the 
fact that not only the pressure drop between the outer and inner boundaries of the spherical shell, but the pressure values 
at these boundaries, can vary within a limited range. In this formulation of the problem, when there is only thermal 
action, the hollow sphere does not completely pass into the plastic state. The research results provide predicting the 
behavior of an object (a hollow sphere) that experiences centrally symmetric distributed power and thermal external 


influences. 


Keywords: hollow sphere, thick-walled spherical shell, thermoelastoplastic state, equivalent stress, associated plastic 


deformation law, stress hodograph, model behavior control parameters. 
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Introduction. The solution to the problem of a thick-walled spherical shell experiencing different external 
influences is given in the monographs [1, 2] and a number of scientific papers on the theory of elasticity, plasticity, and 
thermoelastic plasticity [3-9]. Usually, the case is considered when the loading process is simple. 

The problem of a thick-walled spherical shell is one of the simplest elastoplastic problems when the fields of 
external actions and internal parameters have central symmetry. Due to the central symmetry in the plastic region, the 
regime of complete plasticity 1s performed. For an ideal plastic body, the problem is statically definable, which allows it 
to be solved under any plasticity conditions. In the monograph [1], the most complete solution to the spherical shell 
problem is given, when the dependence of the plastic limit on temperature is not taken into account, and the plasticity 
condition does not depend on the first invariant of the stress tensor and the sign of the third invariant of the stress 
deviator. The cases of only thermal and combined loading are considered, when the temperature on the walls of the 
sphere is set, the pressure on the inner walls is set, and there is no pressure on the outer wall. In [10-13], 
thermoelastoplastic state of various objects was studied, and in [3-9], the process of thermal loading and unloading of a 
sphere free from external forces and a hollow sphere for the Tresca condition was considered with account for the 
dependence of the plastic limit on temperature. The solution to this and similar problems is of interest since it is 
possible to obtain an analytical or partially analytical solution for various mathematical models. An analytical solution 
can be obtained through selecting piecewise linear plasticity functions [11, 13]. Mathematical modeling of objects 
enables to predict their state and behavior depending on the values of the initial parameters [14, 15]. 

Materials and Methods. Problem Statement. We consider the problem of a thick-walled spherical shell (a 


hollow sphere) experiencing centrally symmetric external influences: pressure p, on the outer wall at p=b5 and 
pressure p, on the inner wall at p=a. The thermal effect on the sphere is also considered: temperature 7 is 
maintained at the boundary p=a, temperature 7) is maintained at the boundary p=D. It is assumed that the sphere 


exhibits elastic and plastic properties. The desired state parameters at each point of the sphere are the components of the 
stress tensor, the components of the strain tensors, and the displacement vectors. In the elastic state region, the elastic 
deformations are complete (there are no residual deformations). 

Basic Ratios. All relations are reduced to a dimensionless form. The outer radius of the sphere b is selected as a 
length scale. All values having the stress dimension are assigned to the plastic limit under uniaxial tension k . The scale 
unit for temperature is | C. 


Due to the specified symmetry of external actions, in the spherical coordinate system p,0, of the matrix, the 


components of the stress and strain tensor will have the form: 


o, 0 0 é, O 0 
(s)=| 0 o, OF}, (e)=| 0 &s O 
0 0 Oo, 0 0 6, 


In this case, the equalities 6, =0,, €&, =, are fulfilled. 


o 9 
If the plasticity functions do not depend on the first invariant of the stress tensor and the sign of the third 
invariant of the stress deviator, then, when solving the sphere problem, the plasticity functions will be reduced to the 


form: 
f =|9,-9, Ek. (1) 


Consider the plasticity condition: 


e((G,-0e,) 4 200,06.) 2 


0.,0,, &,&) J= 
fl p ts) p ay) c+n(l+a)” 


(2) 
_ Mo, 5, — HE, —& JI" FAG, — 5 — OE, =6,))" 


ct+n(1+a)” 
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When the parameters have the values: ¢=0, n=1,6=0,a =0, m=1, u=0,k =k,, the condition (2) implies the 
condition (1). In Fig. 1, the plane o,,0, shows the plasticity curves determined from the formula (2) for different 


values of the numerical coefficients in the plasticity function. 





a) b) 
Fig. 1. Plasticity curves: a) for parameters: ¢ =0,2; w= 2; 1 =0,5; 6 =0; m =3; k =1; solid line a =0,5; 


dotted line a@ = 0; 5) for parameters: o =0;¢=0;5=0;4 =1 


The results presented in Fig. 1 show that when the first invariant of the stress tensor is taken into account, the 
radial and circumferential stresses can vary in a limited range when the point of the sphere is in an elastic state. 
Accordingly, the pressure on the boundaries of the sphere should also be limited. When the first invariant in the 
plasticity condition is not taken into account, the elastic state is possible for any pressure value at the boundaries of the 
sphere, but the pressure drop is limited Ap = p, — p,. Taking into account the sign of the third invariant of the stress 


deviator, as noted above, affects the values of the plasticity limits. 
If the values of the state parameters 0,0, determine the point of the region bounded by the plasticity curve, it is 


assumed that the defining equations connecting stresses and deformations are the relations of the Duhamel-Neumann 
law [1, 2]: 
Ee, =(l-v)o,-vo,+ kal, Ee, =0,—-2vo,+E£a?, (3) 
where the Young's modulus F' and the Poisson's ratio v are constants. 
If the state parameters 6,,0, determine the points on the plasticity curve, then an additive representation of the 
total deformations in terms of reversible and irreversible deformations is assumed: 
E,=E, +, € =E +E. (4) 


Complete deformations are determined through displacements from the formulas: 





u du 
2) oo (5) 
p dp 
Complete deformations are bound by the condition of compatibility of deformations: 
d 
ig “8 4+¢,-8 =0. (6) 
r 


Increments of irreversible deformations are related to stresses by the normal law: 
ds’ 
de 2, (7) 
Of/0c, of /0o, 
The relation (7) is generally non-integrable when choosing nonlinear plasticity functions [16]. In the quasi-static 
approximation, the stresses must satisfy the equilibrium equation: 





do, 
fe Ip +2(6,—-0,)=0. (8) 


Equivalent Stress. The equivalent stress is the convex isotropic scalar functions of the stress tensor. In special 
cases, the term “equivalent stress” is synonymous with other terms, for example, “stress intensity” [17]. In this paper, 
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the equivalent stress coincides with the plasticity function. In this case, the equivalent stress will not have a 
discontinuity at the elastic-plastic boundary. 
Temperature Field. The temperature field in the sphere is found from the solution to the boundary value 
problem [1]: 
dT dT 
pt 0, 
dp dp (9) 
T la = ie T p= fie 





The solution to the problem (9) is presented in the form: 
Pe a re. (10) 
(b—a)\p 
Elastic Area. In the region of the elastic state of a hollow sphere, the formulas for stresses have the form: 
BAX B KX abEaAT 
0-Ar -—, 0, =A-_2-_, ha 
pp 2p 2p (l1—v)(b—a) 
Plastic Area. We select the conditions (1). Let us consider the case of only the thermal effect (10). Then, the 
plastic region will originate at the inner boundary of the shell under the condition [1]: 
Ea 2(a° +ab+b’)k 
p=" jarpp =A ter 
(l—v) b(a+2b) 


Denote by c, — the radius of the elastic-plastic boundary p =c,. During the loading process, when B > B,, the 











plastic area a<p<c, increases. When the condition (1) is selected, the stresses in the region a<p<c, are calculated 
from the formulas: 

o! =2KkIn(p/a), 6,’ =o +Kk, K, =sign(o,-9,), 
where K, = sign(o, —9,)|,_. - If AT >0, then «, =—1, if AT <0, then «, =+1. 

If the region c, <p <b remains elastic, then the values A, B and radius of the elastic-plastic boundary c, are 
determined from the conditions of continuity of stresses at the elastoplastic boundary, and the boundary condition 
o|,=0. So, if A and B are determined only from the conditions of continuity of stresses at the elastoplastic 
boundary, then, the following expressions take place: 


2: 
At Oe, Ba ee. (11) 
a) 3) 3¢, = 2 





The equation to calculate c, will have the form: 


3 2 
2,k) In 4). 3-26 fp pare) eer (12) 
a} 3 3b 3c, 3b Ob 


If A and B are determined from the conditions of continuity of stresses at the elastoplastic boundary, and the 








conditions o|,,=0, then the following expressions take place: 
2 3 2 
fi) es eee Ge (13) 
3b° )b ~~ 3b 3 3 
The choice of formulas (11) or (13) affects the steps of the algorithm for solving the problem, but does not affect the 
final results. 








A second plastic region will be generated at the boundary p = b if the following condition is met: 
(0, —6,) L2= kk, K, =k. (14) 
To determine the value AT = AT, when the condition (14) is satisfied, it is required to combine the system of equations 
(12), (14). Since the parameter B enters the equations (13) and (1) linearly, it is possible to obtain a separate equation 


for determining the radius of the elastoplastic boundary: 
2xykin| S.) +O aNG —P) —(0 


9 


a 3(b+¢, be, (15) 
as well as the formula for calculating the parameter Bp : 
2k, («,b° —«,c, )(b—a) 
B =p, = — (16) 


(b° —c’)ab 
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Under further loading, when the inequality B >, is satisfied, the sphere region c, <p <b goes into a plastic 


state at the boundary p=). 
A hollow sphere when exposed to heat. Consideration of material hardening. Consider the case when the 
plasticity condition has the form: 
(6,5), €&,€) )=|6, —9, —d(e; —Ef J =k. (17) 
If there are no residual deformations in the sphere before loading, then as a result of thermal heating, the plastic 
zone will be generated at the inner boundary p=a, when the conditions (15), (16) are met. With further loading, 


plastic region a<p<c, is formed. To find the stresses in this region, it is required to get the corresponding equations 


from the system of equations (3-6), (8), (17): 














2 d°o, do, 6 dabEaAT 
0 5 + 4p — —__——_| k + —————_ |= 0, 
dp dp 1+3d0-v) o(b-a) (18) 
p do, 
Oo, =— O.. 
2 dp 
The solution to the system (18) is written as: 
6, = [thin teen |S 7 
+ —v a 
‘ (19) 
_ 1 Lop _ 36abEaAT ree 
1+36(1-v) -2p(b—a) a) 


The values C,, C,, included in the formulas (19), are determined from the boundary condition o, |_.=0 and 
the condition for the absence of plastic deformations at the elastic-plastic boundary p=c;: 


C, =2«,(k-k,)c -3EaN,c, , 


2 — E : 
C, =2«,k, ing AIG 5 FON 32) (20) 


3a a a 


7 dabAT 
> +38(1-v))(b-a) 
As a result of substituting (20) in (19), we get: 


— 1 retin! 4 — ae 
* -14+36(1-v) a 


_( 2K,ke; fi dabc; EaAT \f 1 1 ri 2«, ke, _ 1 
3 b-a ap 3 la py 

Oo, = a Kk + ee ns 
1+306(1-v) 


a b-a 
Kkce’ dabc EaAT \( 2 1 Khoo ( 2 1 
- +$ ——+ ——— || —+— | |+—| —+— |. 
3 2(b-a) ap 3 la p 


From the solution to the elastic problem, it follows: if AT >0, then «x, =—1; if A7’ <0, then «, =+1. 





a 20 





Accounting for plastic compressibility. Consider the case when the plasticity function is linear with respect 
to the components of the stress tensor: 
G(20, +9, )+K(o, —9,)+N(0, —9,) =k, 1-x7), 


Kk = sign(O, —9O, ). 


(21) 


The plasticity condition (21) can be represented as: 
ao, +Bo, =kU-xT7), 
a=2¢6+K-N, 
B=¢-K+M,. 
Taking into account the introduced notation, to determine the stresses in the plastic region, we obtain the 


problem: 
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ao, +Bo, =kI- x7), 





do, 
p+ (9, -9,)=0 (22) 
6, | a= DP, 


The solution to the problem (22) has the form: 


: (M-l)k, 2Nk, (2) ~ -1 2N | 
o =|—p, +—*+—_ "+ + —__ ** ]}|- —| —_+_+——_4+_|f,, 
: ° a+p a(a+2B) /\p a+B p(a+t+2B) 


iB (M,-l)k, 2Nk, \(a\y*° (M,-1 WN, 
o, =-—| —p, +—*>*—_ + —*—__ | - — + ——*— |k,, 
(oF a+ a(a+ 2B) )\ p a+B p(a+2p) 
aAT N= abAT 
b-a’ b-a- 


To get the correct result from (20), when, for example, a+f=0, it is required to perform a limit transition 


(23) 











where the notation is introduced: M =T, — M, = YM, N=. 


when solving (22). It is easier to obtain the correct result directly in (23), while taking into account that the condition 
a+ =0 is met. In this case ¢ =0, so we get: 


2 
Go, = Ke ‘aan [4-1], |p , 
Oo a\p a 


ua d- )0-in®) +, (2-2) 9, 
a p 


as 
Ou a 





During the loading process, the plastic zone originates at the boundary p =a, when the following condition is 
met: 
AT =kt=1,)| bar SSPE ONES) 
2(1—v)(a +ab+b ) 
Research Results. Fig. 2 shows stress graphs and stress vector hodographs, when the sphere region 
corresponding to the condition a<p<c,, is in a plastic state, and the sphere region corresponding to the condition 


c, <p <b, is in an elastic state. 


Oo 


» ea Oo, =—,e€a 








M45 
a) b) 


Fig. 2. Stress graphs (a) and stress vector hodographs (5) 
for parameter values:: k = 1; v = 0.3; a=0.5; b=1; AT =170; c, = 0.57 


Fig. 3 shows stress graphs and stress vector hodographs, when the sphere regions corresponding to the 
conditions a<p<c, and c,<p<)b, are in a plastic state, and the sphere region corresponding to the condition 


c, <p <c,, 1s in an elastic state. 
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fo) 
Oo, =—.ea 
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a) b) 


Fig. 3. Stress graphs (a) and stress vector hodographs (5) for parameter values: 
v = 0.3; a = 0.5; AT = 270; c, = 0.62; c, = 0.88 


Fig. 4 shows stress graphs and stress vector hodographs, when the sphere region corresponding to the condition 


a<p<c,,is inaplastic state, and the sphere region corresponding to the condition c, <p <5, is in an elastic state. 


fe) 
—,ea 





A5 


a) b) 


Fig. 4. Stress graphs (a) and stress vector hodographs (4) for parameter values: 
v = 0.3; a =0.5; AT = 215; c, = 0.58 


Fig. 5 shows stress graphs and stress vector hodographs, when the sphere regions corresponding to the 


conditions a<p<c, and c,<p<)b, are in a plastic state, and the sphere region corresponding to the condition 


c, <p<c,, 1s man elastic state. 
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a) b) 
Fig. 5. Stress graphs (a) and stress vector hodographs (5) for parameter values: 


v = 0.3; a=0.5; AT = 270; c, = 0.61; c, = 0.86 
Fig. 6 shows stress graphs and stress vector hodographs when the sphere regions corresponding to the 
conditions a<p<c, and c,<p<)b, are in a plastic state, and the sphere region corresponding to the condition 
c, <p <c,, is in an elastic state. 


,ea 


. 
1 





a) b) 
Fig. 6. Stress graphs (a) and stress vector hodographs (6) for parameter values: 


v =0.3; a=0.5; AT = 79; ¢ = 0.1; n = 0.1; Ea = 0.012; y = 0.0017 


Discussion and Conclusions. The calculation results show that in this formulation of the problem, when there 
is only thermal action, the hollow sphere does not completely go into the plastic state (Fig. 2-6). Hardening causes an 
increase in the equivalent stress in the plastic region and a decrease in the radius of the elastoplastic boundary (Fig. 4, 5). 
The elastic region cannot completely disappear under loading. Plastic compressibility and the dependence of the plastic 
limit on temperature have a significant effect on the stress state of the hollow sphere (Fig. 6). 
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Introduction. A review of the application of theoretical mechanics methods for the development of algorithms of 
approximate analytical simulation of a motor vehicle (MV) movement in case of oblique collision with side cable 
barriers is performed. The representation of the MV as a material point in this type of collision is validated. The 
study objective is to demonstrate the application of a physical model of a material point to describe the motor vehicle 
dynamics in the event of its oblique collision with side cable barrier. 

Materials and Methods. A new physical model that describes the opposition to the motor vehicle movement from the 
side of a cable barrier in an oblique collision is proposed. New methods of approximate analytical construction of the 
MV movements during an oblique collision with the side cable barriers are presented. The analytical calculation 
results are verified by the data of the finite element (FE) simulation of the collision according to the data of field 
tests. The FE simulation was carried out using a multi-purpose finite element complex LS-Dyna. 

Results. New analytical algorithms have been developed for the MV movement in case of an oblique collision with 
side cable barriers, as well as a new physical model describing the opposition to the MV movement from the side of 
cable barriers. The application of a physical model of a material point to study the motor vehicle dynamics during an 
oblique collision with side cable barriers is established scientifically, including the comparative analysis of the 
kinematic results of the virtual test with kinematic calculations obtained on the basis of algorithms for analytical 
construction of the MV movements. 

Discussion and Conclusions. The analysis of the kinematic results of the virtual test in comparison with the 
analytical kinematic calculations has shown that the representation of a motor vehicle as a material point in case of 
an oblique collision with side cable barriers 1s reasonable since the MV movement is close to translational motion. 
The results obtained can be used in the development and analysis of the correctness of the FE modeling of a side 


collision of a motor vehicle with cable barriers. 


Keywords: cable barrier, analytical modeling, motion trajectory, displacements, material point, motor vehicle, oblique 


collision. 
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Introduction. Safety systems designed for the organization of traffic flows on road routes help to minimize the 
adverse consequences of road traffic accidents (RTA). Currently, cable road barriers are gaining popularity. They are 


installed on the road dividing strip to prevent intentional and unintentional crossing of the road by vehicles, as well as 


Kolesnikova G. P. Material point physical model rationale while studying kinematic characteristics of a motor vehicle 





on the roadsides, to prevent the departure of vehicles outside the road’. Different structures of cable barriers consist of 
common elements: cables, racks, anchor blocks, clamping devices, but differ in the type of their installation’ [1-4], in 
particular, in the method of fastening the cables (Fig. 1), which significantly affects the development of mathematical 
models. The behavior of the racks under the MV collision is also different. Racks can be crushed or, coming out of the 
sleeves, assembled with one another, and thus provide additional braking of the MV [5]. 

The collision of a car with a cable barrier is characterized by complex interaction mechanics, since the cables 
have a high degree of geometric nonlinearity, and the racks and the ground have a high degree of physical nonlinearity, 


while almost all processes are transient [5]. 





Fig. 1. Some types of cable racks [6] 


To determine the MV trajectory, string vibration equations can be used in the analytical calculations of 
arresters’ * and cable-stayed structures” °. But due to the complexity of the construction of barriers, there are difficulties 
in the formulation of boundary conditions. Energy methods also give a set of equations that are not solved analytically’ 
[7]. All this causes the need to study barriers with the help of engineering software packages for both modeling the 
cable barrier itself and the crash test system. In this case, the numerical finite element method in the explicit formulation 
is often used, which is implemented in the software complexes LS-DYNA, MARC, NASTRAN, etc. [5, 8, 9]. 

The disadvantages of building FE models are the following: the development complexity; the duration of the 
calculation period (120 hours or more); the need to check the correctness of the construction. Thus, the period of 
preparation of the basis for the study is quite long. In this regard, a method of approximate analytical calculation of the 
MV movement from the data of field tests, based on the methods of theoretical mechanics, is proposed. The analytical 
calculation provides checking the correctness of the developed FE model, as well as reducing the time of the passive 
safety study by 2-3 times, since for the analysis of passive safety, it is no longer necessary to build a FE model of the 
cable barrier itself, but to use the analytically obtained movements. 

Materials and Methods. As a result of observing the MV during full-scale crash tests, a hypothesis was put 
forward: in case of the collision with the side cable barriers, the MV can be represented as a solid body making a 


translational motion (Fig. 2). 


' Industry road methodological document ODM 218.6.017-2015. Guidelines for the use of road barriers of various types on federal highways, 
recommended by the order of the Federal Road Agency, dated December 23, 2015, no. 2489-r. URL: 
https://files.stroyinf.ru/Data2/1/4293757/4293757596.pdf (accessed: 08.01.2021). (In Russ.) 

* Home DA. Report 350 Acceptance of New York Three-Strand Cable Terminal. Office of Highway Safety Infrastructure, FHWA, U.S. Department 
of Transportation — Washington, D.C.; February 14, 2000. Available from: https://highways.dot.gov/ (accessed: 21.01.2019). 

* Nuralieva AB. On dynamics of the space elevator cable: Cand.Sci. (Phys.-Math.), diss., author’s abstract. Moscow, 2012. 20 p. (In Russ.) 

* Mikhailyuk DS. Finite element modeling and research on the dynamics of a deck arrester: Cand.Sci. (Eng.), diss., author’s abstract. St. Petersburg, 
2009. 19 p. (In Russ.) 

> Dyadkin SN. Rationale, technology of balanced erection and monitoring of byte bridge spans with account for climatic factors (case study: the 
bridge over the Ob River near Surgut): Cand.Sci. (Eng.), diss., author’s abstract. Volgograd, 2005. 20 p. (In Russ.) 

° Le Thu Huong. Optimization of parameters of spans of suspension bridges, reinforced and not reinforced with inclined cables, when designing them 
with using PC: Cand.Sci. (Eng.), diss., author’s abstract. Moscow, 1999. 25 p. (In Russ.). 

™ Los MV. Numerical modeling of the behavior of the “body-rope” system with account for flexural rigidity of the cable and mechanism of looping: 
Cand.Sci. (Phys.-Math.), diss., author’s abstract. Moscow, 2000. 19 p. (In Russ.) 
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trajectory of MV points 


Fig. 2. The model of MV as a solid body performing translational motion 


In addition, a number of fundamental assumptions on the dynamics of the MV during an accident are made, 
justifying the use of the proposed theoretical and mechanical models of the phenomenon: 

— the lateral and longitudinal decelerations of the vehicle are constant during the time interval required for the 
MV to be oriented parallel to the undeformed barrier; 

— vertical and rotational accelerations of the vehicle are ignored; 

— the lateral component of the speed is zero after the vehicle is redirected parallel to the guardrail; 

— as the vehicle is being redirected, it does not engage with the guardrail; 

— the deformation of the vehicle occurs in the collision zone, but its center-of-gravity position does not 

change significantly; 

— the MV center of mass moves as if all its mass is concentrated at this point; 

— the barrier can be rigid or flexible; 

— the friction forces of the car tires on the road surface are ignored; 

— the guardrail system does not contain breaks that can cause sudden vertical movements of the vehicle. 

With this representation of motion, a MV (Fig. 3) can be considered as a material point of some mass m [10]. 
The Cartesian coordinate system origin corresponds to the point of origin of contact between the vehicle and the barrier. 


The motion of the material point begins with at the speed Vo, directed at an angle a to the plane of the fence (x-axis) [6]. 


y 





Fig. 3. Calculation scheme of MV [6] 
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The opposition to the movement of the vehicle from the side of the cable barrier was modeled by the 
forces [6]: 
— with lateral deviation along the y-axis — by the elastic force P =—Ay and dissipative resistance 
P* =—by ,ie., Fy, =—ky — by; 
— the opposition to movement along the x-axis — by the force of friction / =—fN,, where N, — the force 
of inertial pressure, 1.e., NV, =my [6]. 
As a result, the MV motion is described by a system of differential equations [6, 11]: 
t= — fy; 
yr+2ey+ py =0. 
The solution to the system is displacement as a function of time f: 


x(t)=—fy(t)+v, (cosa+ f sina); 


y(t)= yee sin pt. 
P 





In [6], it is proved that the Coulomb friction force from the side of the road barrier does not significantly affect 
the simulated the MV motion. Physical modeling of the MV as a material point is not new under studying the case of 
collision with lateral barriers. This approach was used by R. M. Olson, E. R. Post and F. F. McFarland when describing 
a car hitting a rigid bridge guardrail [12]. Here, the calculation of the side barrier resistance was based only on the 
classical Amontons-Coulomb model of friction. 

In contrast to the problem, the nature of the interaction between the MV and the barrier is described in this 
paper by a fundamentally new model. With the external similarity of the problem statement, the physical essence of the 
interaction between the MV and the barrier differs qualitatively from the essence of the interaction with more rigid 
fences. In our paper, this was taken into account through the Amontons-Coulomb force models considering the force of 
the inertial normal pressure and the Kelvin-Voigt resistance model. 

Also, in the context of the study on the nature of the impact on cable barriers, the paper by M. B. Bateman and 
others should be mentioned [9]. The data of full-scale tests presented in this work demonstrate clearly that in case of 
oblique collision with cable barriers, the MV motion is close to translational one, when the yaw angle does not exceed 
10°. Here, the head-on crash process is described by two models: 

1. A simple dynamic model of a vehicle, where the Runge-Kutta fourth-order method was used for the 
numerical solution to the differential equation of the MV motion. In this case, the forces of action from the side of the 
barrier (cable tension) and the road (friction force according to the classical Amontons-Coulomb model) are considered. 
The wheels are not assumed to require steering under impact. 

2. A quasistatic model of the barrier, which is designed to calculate the change in tension of the cables when 
their geometry changes as a result of the movement of the car. 

As a result, the computational model is quite cumbersome and requires a numerical solution, whereas the 
model proposed in this paper provides calculation algorithms that are easily carried out analytically without the use of 
numerical methods. 

Approximate analytical calculation of the MV kinematic characteristics. As a result of the assumptions 
made, an algorithm for constructing the MV motion called harmonic has been obtained [6]. Further, we propose 
quadratic and cubic algorithms based on the construction of polynomial displacement functions [13]. 

Here, the opposition to the MV motion from the side of cable barriers along the y-axis is determined by some 
function P,(t), the result of its integrating is as follows: the MV displacement function along the y-axis has the 
character of a second or third order polynomial. The simulation of the opposition to the movement along the x-axis 
remains the same. 

We will conduct a comparative analysis of the results of the approximate analytical calculation of the MV 
kinematic characteristics. Fig. 4-6 show the computation data for bus Mercedes-Benz 0345. When constructing the 
movements of the MV, data of field tests of the State Research Center of the Russian Federation FSUE “NAMI” were 
used. 
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quadratic approximation — dotted line; harmonic approximation — gray line [13] 
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Fig. 6. Dependences of the deviations x, m, and y, m, on time ¢, s, calculated by the methods of cubic and harmonic approximation 
along the x-axis — black solid line; along the y-axis: dotted line — by the method of cubic approximation, gray line — by the 
harmonic approximation method [13] 


As can be seen from the above calculated data, the divergence of analytical methods does not exceed 20 %. 
The same calculations were carried out for car GAZ-3102 [12]. 

Construction of a finite element model of the contact of the cable barrier and bus Mercedes-Benz-0345. 
To analyze the performance of the approximate method of constructing the MV trajectory, a FE-model of the 
interaction of the cable barrier and bus Mercedes-Benz-0345 was developed in accordance with the industry standard 
STO 521000-005-10690827-2015', agreed with Rosavtodor in 2017. According to this standard, the installation of 
cable barriers for tests was carried out; the test results were also used in the analytical modeling of the trajectories and 
deviations of the MV. The object of virtual tests were cable barriers of 14DD/U 4(300)-P-1, 1-3, 0-GB brand [11, 13] 
with the following parameters: 

— octagonal sleeves of GZ 500/U brand; 

—STD-2 racks consist of two elements in the form of a square pipe with a length of 1500 m with a cross- 
section size of 50=50=3 mm; 

— spacing of racks with concreted sleeves at the working area — 3 m; 

— spacing in the initial and final sections — 2 m; 

— the height of the racks above the road surface — 1.1 m. 

Cable parameters: 

— three-strand, seven-conductor; 

— diameter —19 mm; 

— rated breaking force ~ 173 kN; 

— the number of branches — 4. 

Bus Mercedes-Benz-0345 parameters [12, 13]: 

— gross mass — 14,050 kg; 

— overall length — 12,000 mm; 

— overall width — 2,500 mm. 


The speed of the MV coming into contact with the barrier 1s 69 km/h. The angle between the plane of the 
barrier and the direction of the MV movement — 20° in accordance with the methodology of GOST R 52721-2007 


“STO 521000-005-10690827—2015. Retaining side barriers for vehicles with the use of a C-section of a beam. “PIK” Enterprise. rosavtodor.gov.ru 
URL: https://rosavtodor.gov.ru/storage/app/media/rosavtodor/b/2016/02/09/sto_ 521000 005 10690827 2015.pdf (accessed: 10.01.2021). (In Russ.) 
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(nos. 6.1-6.6) [11, 13]. Time-lapse images of the MV according to the data of full-scale’ and virtual tests are presented 
in the Table 1. 


Table 1 
Comparative frame-by-frame images of MV 
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*STO 10690827-001—2015. Lateral deformable road barriers, cable type. Specifications. “PIK” Enterprise. rosavtodor.ruwwwrosavtodor.ru URL: 
http://rosavtodor.ruwww.rosavtodor.ru/storage/app/media/uploaded-files/sto-10690827-001-2015.pdf (accessed: 08.01.2021). (In Russ.) 
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The field test results practically coincided with the results of the virtual tests, since in both cases, the dynamic 
deflection of the barrier was 1.5 m; the path length of the interaction of the MV and the barrier from the calculation was 
23.5 m, according to the field tests — 25 m. The relative error of the path length of the interaction did not exceed 6 %. 

Research Results. As a result of a virtual test using the developed FE-model of the cable barrier, the 
trajectories of the characteristic points of bus Mercedes-Benz-0345 were determined. Fig. 7 shows the points of the MV 


for which the readings were taken, point C corresponds to the center of gravity. 





Fig. 7. MV points selected for analysis: A, B — on the frontal surface; 
D, E — on the rear surface 
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Fig. 8. Dependence of displacements of the MV characteristic points on time. 
Dotted line — for points A and B along the y-axis; solid lines — along the x-axis: 
black — for point A, gray — for point B 
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Fig. 9. Dependence of displacements of the MV characteristic points on time. 
Dotted line — for points D and E£ along the y-axis; solid lines — along the x-axis: 
black — for point D, gray — for point EF 
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As can be seen from the presented data, the displacements of points A, B, and D, E along the y-axis match in 
pairs. The same can be observed for speeds. Thus, in case of oblique collision of the MV with the side cable barriers, 
the movement of the MV is really close to translational. 

Discussion and Conclusions. The analysis of the kinematic results of the virtual test in comparison with the 
kinematic analytical calculations shows that the representation of the MV as a material point in case of its oblique 
collision with the side cable barriers is justified, since the movement of the MV is close to the translational one. 
However, if necessary, for a more accurate analytical calculation of the kinematic characteristics of the motion of points 
located in areas close to the points C, D, FE, an additional analytical calculation is required, for example, using the 


equations of plane-parallel motion of a solid body. 
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Introduction. When organizing drilling operations, one of the major problems is the accuracy and smoothness of 
lowering bundles of pipes into the shaft of the drilling rig. This depends on many factors, including the operation of the 
hydraulic brake of the lifting device. The objectives of this work are to create and study a mathematical model of 
hydrodynamic ballast in a drilling rig. Using the inverse analysis method, the effect of some performance indicators on 
the braking torque of the hydraulic brake is studied. 

Materials and Methods. The experiments were performed using a laboratory setup, which is a model of a hydrobrake. 
Its valve was closed under various conditions to obtain several pressure values with the calculation of the braking 
torque when a certain weight was suspended. The real (field) operating conditions of the hydromatic brake were 
simulated, and the results obtained were compared. When creating a mathematical model, the inverse analysis method is 
used. It is based on the results of experimental measurements and provides expressing the totality of the effects of 
individual variables on the braking torque. 

Results. A mathematical model of the hydraulic brake has been created and tested. The dependence of the braking 
torque on the pressure, density, and viscosity of the ballast fluid is determined. The influence of each variable is 
determined experimentally since the dependence under consideration cannot be represented as a direct relationship. The 
inverse analysis method is used to obtain a set of constant values that give the optimal solution. Taking into account the 
standard error array and the minimum standard error, the statistical errors made during experimental measurements are 
considered. The physically acceptable range of values of the proposed mathematical model is visualized. Using a basic 
(nonlinear) mathematical model, the auxiliary braking torque of a hydrobrake is calculated as a function of pressure, 
density, and viscosity. The proposed model validity is established. The calculated values of the braking torque were 
used as a criterion of correctness. The erroneous discrepancy did not exceed 6 %. For additional testing of the model, a 
computational experiment simulating field conditions was performed. 

Discussion and Conclusions. For mathematical modeling of hydrodynamic ballast in a drilling rig, it is advisable to use 
the inverse analysis method. The model proposed in this paper relates the braking torque of a hydrobrake to the 
Operating parameters of the fluid inside the ballast: pressure, viscosity, and density. The objectivity of the model is 
validated. An amendment to it is proposed to simulate the operation of the brake in the field. Based on the results 


obtained, in future studies it is advisable to test the created model in the field with a real payload. 
Keywords: hydromechanical ballast, mathematical modeling, inverse analysis method. 
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Introduction. When setting up drilling rigs, the operations of lifting and lowering drill pipes are provided with 
lifting gears (coils). They are equipped with two types of brakes, the main of which is mechanical. During the drilling 
process, the weight of the pipes exceeds the payload. To compensate for it, an auxiliary control ballast designed to 
reduce the speed when installing a group of drill pipes is introduced [1-4]. 

Ballasts are used to control the load and speed of the hook, as well as to absorb the kinetic energy of a group of 
borehole drill pipes. In addition, ballasts: 

— reduce the effort on the drilling rig, especially under heavy loads; 

— reduce the wear of the main mechanical ballast elements; 

— help to slowly and smoothly stop the load attached to the hook. 

One of the types of auxiliary brakes is hydrodynamic. In this case, the water in the ballast converts part of the 
absorbed kinetic energy of the lifting axis into heat, and due to this, the pipes are lowered and raised [5—7]. 

When manipulating the pipes, the brake axis is connected to the axis of the lifting gears. During the operation, 
the moving part of the installation displaces water in the direction of the inclined blades inside the stator, and rotates at 
a speed equal to the rotation speed of the axis of the lifting gears. Water falls on the stator ribs, and then on its parts. 
The process is repeated, and the forces that hinder the movement of the rotor increase. This creates a braking torque that 
reflects the rotational movement of the lifting gear axes and reduces the rotary speed. As a result, speed of lowering the 
group of drill pipes decreases [8—10]. 

The braking safety factor (if its value does not exceed the permissible value) coordinates the auxiliary 
hydrodynamic work and the effects of the operation of two main brake systems, as well as provides a longer service 
time for the main brake elements. This 1s done through reducing the wear of the friction discs and the planes of the 
brake wheels. Increasing the braking torque of the ballast provides the correct braking movement of the cylindrical 
hoist. This demonstrates the importance of studying the hydrodynamic brake. 

The hydraulic power of the brake N is determined from the ratio: 

N=pig*O-i, (1) 
where p — the density of the working fluid (fresh water); g — gravity acceleration; OQ — the amount of working fluid 
consumed, equal to the volume of the working fluid, which passes through the system of ribs in one working cycle; H 
— the height (level) of the working fluid in the ballast. 

The braking torque Mz of the forced engine braking is determined from the ratio: 

M,=p:g:O-H/o, (2) 
where m — angular velocity of rotation of the moving section of the ballast. 


Indicators of kinematic braking of ribbed hydraulic machines are determined from the known theoretical ratios 


of the braking torque: 
M, =i,,-p:(D’-d’)-o’, (3) 
n- 
M, =i, -p-(D° —d’ )-—. 4 
aay “P-( ) a0 (4) 


Here, D — the outer diameter of the “ring” of the working fluid formed during the rotation of the rotor wheel and 
assumed to be equal to the diameter of the wheel of the moving section; d — the inner diameter of the working fluid 


ring, which depends on the level of the ballast; » — the number of revolutions of the moving section of the ballast (the 


2mn > on 
o=——>@ = 
60 100 


number of rotating axes), rpm: 





Antypas I. R., et al. Inverse analysis method for mathematical modeling of hydrodynamic ballast in a drilling rig 


The hydraulic braking torque coefficient 1,, 1s a dimensionless value that takes into account the shape of the 
working cavities in the ballast, the parameters, and the number of brake ribs. In practical calculations, the average value 
of the hydraulic braking torque coefficient is assumed to be 0.3. 

The objective of the study is to investigate the effect of some performance indicators on the braking torque of 
the hydraulic brake, which is the pressure inside the brake chamber, as well as on the density and viscosity of the 
working fluid inside the ballast. These indicators are absent in the ratio (4); therefore, we assume that the brake operates 
at an atmospheric pressure of 1 atmosphere, the fluid is fresh water with a density of 1 g/cm’ and a viscosity of 1 stoke. 

Materials and Methods. Laboratory experiments were conducted at the University of Aleppo (Syrian Arab 
Republic) on a device that is a model of hydrodynamic brake inhibition (Fig.1). 














Fig. 1. Laboratory installation — hydrodynamic brake inhibition model: 

1 — pulley, 2 — pulley armrest, 3 — pressure gauge, 4 — exhaust valve, 5 — thermometer, 6 — opening in the upper part of the 
tank, 7 — fluid level control valves, 8 — fluid outlet channel, 9 — fluid outlet from the tank, 10 — fluid drain valve, 11 — inlet 
valve, 12 — inlet line, 13 — hydraulic brake, 14 — coil, 15 — cycle number meter, 16 — cable, 17 — payload, 

18 — communication hub, 19 — control line, 20 — fluid tank. 


Stages of experiments 

1. The tank and ballast are filled with fluid to the required level (0.106 m to the level of the first valve). 

2. Using a hand lever connected to the coil, the weight suspended on the hook is raised by 0.317 m. It weighs 8 
kg and is connected to the cable. 

3. The bundle of pipes is allowed to fall under its own weight. 

4. Readings from the indicator of the rotation of the axis of the coil are taken. 

5. Then the impact of the indicators obtained during laboratory experiments on the study of hydraulic braking 
is recorded: 

— the pressure created in the brake chamber (P); 

— the density of the working fluid (p); 

— the viscosity of the working fluid (u). 

First, the braking torque of the auxiliary hydrodynamic brake is calculated from the ratio (4). For this purpose, 
the values n, D, d are determined. 

The outer diameter of the “ring” of the working fluid D is formed when the wheel of the moving section (rotor) 
rotates and is assumed to be equal to its diameter. The diameter D of the propulsor of the laboratory unit in the practical 
experiments is 0.33 m. 

The inner diameter of the working fluid ring d depends on the level of the ballast. To determine this diameter, 
the following volumes are aligned: 

— working fluid inside the ballast; 

— the fluid in the “ring”, which is formed when the movable section rotates inside the ballast and is limited by 
the height of the coolant in the tank (0.106 m). 

At D = 0.33 m, the working volume of the ballast fluid is 2.85 liters. The calculations have shown that the 
inner diameter d = 0.283 m. 

The number of revolutions n is associated with the suspended load. If the brake does not work, this is the 
number of revolutions of the coiler. If it is separated from the brake, then n is taken from the laboratory experiments: 
with a load weight of 8 kg, n = 240 rpm. 
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We calculate the hydraulic braking moment for the weight of 8 kg attached to the hook. During the experiment, 
the tank was filled with fresh water to the level of the first valve (0.106 m) at the following parameters: pressure 
P= 1 atm, fluid density p = 1 g/cm’, and fluid viscosity 1 = 1 St. Therefore, the moment of hydraulic braking: 


2 
fal =355 H-M. 
100 





M,, = 0.3-1000-(0.33° —0.283°) 


Consider the inverse relationship to (2). First of all, we are talking about the inverse relationship between the 
torque and the angular velocity (m). A decrease in the number of revolutions n by 20% with an increase in the fluid 
density means an increase in the torque by 20%. Consider the permissible value of the braking torque 355 N-m. As a 
result, we get the desired value of the braking torque at the following parameters: pressure P= 1 atm, density 
p = 1 g/cm’, and fluid viscosity 1 = 1 St. 

With regard to this work, we note the following. When determining the braking torque of a hydraulic brake, it 
is required to consider the number of revolutions of the coil connected to the brake axis (1.e., during braking): a decrease 
in the winding speed by a certain percentage means the same increase in torque. 

The effect of the pressure inside the ballast on the braking torque of the hydraulic brake. Fluid pressure 
was created inside the hydraulic brake by closing the valve in the fluid outlet line. We conducted a laboratory 
experiment and calculated the braking torque with a suspended load weighing 8 kg. The results are shown in Table 1. 


Table 1 


Change in the braking torque of the hydraulic brake when the pressure inside the ballast changes 


Tettorgue mit, Nm ——~*Y~st0s| toa [103 | 100 | 95 | 90 


Fig. 2 shows the change in the braking torque when the pressure inside the ballast changes. 
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Fig. 2. Change in the braking torque at different pressures inside the ballast 
The effect of the working fluid density on the braking torque. The experiments used chemicals that 
increase the density of water without corroding the elements of laboratory equipment. 
Various concentrations of substances dissolved in water for the production of a ballast fluid of different 
densities are considered. The conditions are the same: the braking torque is investigated with a suspended load of 8 kg 


(Table 2). 
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Table 2 
Dependence of the braking torque on the fluid density, the number of coil turns, and additives in the ballast fluid 


(Number ofeoitums «Y= —~| tf me | YC 


Sodium chloride salts 


—Sameioridesats SSS 
———CC™~*CiR dug SS 


Food sugar 


A very weak change in the viscosity with an increase in the concentration of the sodium chloride salt was 
observed. We can assume that the viscosity of the fluid is approximately equal to 1 St. 
Fig. 2 shows the dependence of the braking torque on the density of the working fluid when using sodium 





chloride salt and food sugar. 
450 
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Fig. 3. Dependence of the braking torque on the density of the working fluid with the introduction of sodium 
chloride salt (a) and food sugar (5) 
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The effect of the working fluid viscosity on the hydraulic system braking. Various combinations of 
concentrations of glycerin and sodium silicate dissolved in water were used to produce ballast fluids that differ in 
viscosity. The conditions are the same: the braking torque is investigated with a suspended load of 8 kg (Table 3). 


Table 3 
Dependence of the braking torque on the viscosity of the fluid and additives 
in the ballast fluid 


Fluid density 
Number of coil tums a 


Glycerin 


Sodium silicate 


Fig. 4 shows the changes in the braking torque with the change in the viscosity of the working fluid when using 





glycerin and sodium silicate. 
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Fig. 4. Dependence of the braking torque on the viscosity of the working fluid 
inside the ballast with the introduction of glycerin (a) and sodium silicate (5) 
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The almost linear dependence obtained in this way should be checked using a mathematical model. 
During the study of the viscosity index, it was found that the density and viscosity parameters do not depend 
much on each other, and the density is not equal to 1 g/cm? (Table 4). 
Table 4 


Fluid density and viscosity ratio 


2 ee A eee Se 





The density influence coefficient did not exceed 6 %. It was calculated as the difference between the densities 
of a viscous fluid and fresh water. For example, if the viscosity was 30 St, the density was 1.025 g/cm’. Hence, the 


density impact factor: 


A= —.100=2.5%. 


Research Results 

Creating a mathematical model of a hydraulic brake. The dependence of the braking torque on pressure, 
density, and viscosity cannot be represented as a direct relationship. The effect of each variable on the braking torque is 
determined experimentally (see Fig. 2—4). 

The inverse analysis method. In the scientific and reference literature, there are no recorded indicators of the 
relationship of the braking torque and pressure, density, and viscosity. Therefore, the inverse analysis method was used 
[4]. It is efficient for creating mathematical models based on experimental measurements. It can be used to show how 


the combination of the above variables affects the braking torque. This can be expressed by the relation: 
Dare), (5) 


where F' — a function that relates the considered phenomenon D to a set of variable values and a set of values for 
constants of the mathematical model of the phenomenon. 
The inverse analysis provides finding a set of values of the assumed constants of the model C through the 


inverse dependence: 


Por (0 76). (6) 
where D,,, — a set of experimental values for the phenomenon under study [4, 5]. 

It is assumed that a direct solution to the relation (6) is impossible. Therefore, an iterative system should be 
used to get a set of constant values that give an optimal solution. The generalization specifies a set of P values for the 
relations: 

D,=F(P,C), (7) 
ID, -D,|<e, (8) 
where D. — a Set of calculated values for the phenomenon under study; ¢ — the required accuracy in accordance with 


the calculation of the phenomenon under study. 
The proposed method takes into account the statistical errors made during the experimental measurements, 
focusing on the standard error array (S,) and the minimum standard error. Thus, according to the statistical Gaussian 


distribution, the dependence of the density of statistical data: 


1 t -l 
(ACg va constenp{ ~3{(D, —D,,) C;'(D. -D, 1} (9) 
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Here, C, — standard error array. The symbol ¢ stands for the matrix: 


Cl ices tO conan 0 
CNW eactenae a (2) eee 0 |,d, (10) 
Denes 0, 555559059995, (72) 


where n — the number of control points. 


D,.. 1s a function of P, so, the problem is related to determining such a value of P, that D. gives the maximum 

value P,. 
S, =(D.-D,,) -C;'(D, -D,). (11) 
The minimum value S, can be obtained for more than one set of P values. Therefore, a correct set of 


parameters is formed using actual physical modeling of the value of these parameters (Fig. 5) [4, 5]. 





Fig. 5. Determining a physically acceptable range of values (green area): 
S'; — physically acceptable mathematical solution, S,, $3; — mathematical solutions 


Thus, we use a set of primary values for the parameters Py, which are derived from the average load of the 
physical field with account for the standard deviation matrix S>. 


When using new Gaussian distribution, we obtain the following statistical intensity parameter: 


f,(P)=P, = constexp| -3[(P-R) -C; (P-2,)]). (12) 
Here, C,, — standard deviation matrix: 
S; (1) peaeeenoeed: 0 eee 0 
Ce) Oraretescemacccen: Ss; (2) heeeaeee 0 |. (13) 
Oaceteteraeien. O) ccpenes teu S; (r) 


The problem expressed in the relation (6) is solved by inverse analysis. In this case, a set of values of P is 
within the limits indicated in (7) and (8). Hence, it is possible to determine a general region of the parameters P, and P, 
from the Gaussian distribution: 


f (P)=P.-P, =const.exp(—S), (14) 


1 


s=-[(D. -D,) +C; (D,-D,)+(P-R)-C; (P-B) |. (15) 


So, for the solution, it is required to find the maximum or minimum values of S of the function f (P). 
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The only way to find the minimum value of S is to use numerical methods, such as Gauss—Newton. This 
approach is based on the transformation of analytical relations into digital iterative ones with account for the error made 


due to ignoring some restrictions in the analytical relations (Fig. 6). 
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the frequency step: i+1 
D iAP :) 
Calculation of brake 
torque values for a set 
of parameters 
Pi. 

i=i+1 


|De- Dale 


Results display 


Fig. 6. Inverse analysis method algorithm 


At the next stage of the work, the auxiliary braking torque of the hydrodynamic brake was calculated using a 


basic (nonlinear) mathematical model. It depends on the pressure, density, and viscosity: 


M,=f(P, p,u)t+C, (16) 
M,=a-p"+b-p” +d-p"+C. (17) 
The mathematical model (17) is a general nonlinear model if the values n,, n», n3 are not equal to one. The set 
of constants to search are a, b, d, and n;, n>, n3. 
After applying the digital iterative method with an accuracy of 0.001, the parameters presented in Table 5 are 
found. 
Table 5 


Parameter values of the proposed model 


[Poamees —[ om [ m [ = | « | 6) | @ 7] c_ 


According to Table 5, the physical phenomenon under study can be represented as the linear model (7), n, 
n3) =]: 
M, =58.5P + 284.6p +1.29u + 10.62 . (18) 


The mathematical model is based on the results of laboratory experiments and implies the homogenization of 
units of measurement in accordance with the value of constants. 

Adequacy of the mathematical model. The first stage of determining the adequacy of the model is a braking 
test at: 

— pressure P = | atm, 

— fluid density p = 1 g/cm’, 

— fluid viscosity np = 1 St. 
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In this case, the braking was equal to the initial braking torque Mz = 355 N-m. This 1s logical, given the 
experimental values. 

The second stage: three laboratory experiments with an 8-kilogram weight suspended on a hook. They are 
briefly described below. 

The first experiment. We took the maximum values of the variables and the braking torque, and then 
performed the calculation using a mathematical model based on the ratio (18). We determined the error rate through 
comparing the experimental and mathematical results. 

The second experiment. We took random values of variables that are relatively far from the experimental 
values and from the calculated values for the braking torque (according to the mathematical model). We calculated the 
error rate through comparing the experimental and mathematical results. 

The third experiment. We took random values for variables that are relatively far from the experimental values 
and from the calculated values for the braking torque (according to the mathematical model). We calculated the 
percentage of errors through comparing the experimental and mathematical results. 

Table 6 shows the results of validating the adequacy of the mathematical model. 

Table 6 
Results of checking the model adequacy 


Braking torque Braking torque 


(calculation) (mathematical model) 





So, the result validates the model adequacy. First, the level of its fallibility is found to be acceptable. Secondly, 
the identified errors have a scientific explanation. The fact is that the density parameter does not depend on the 
viscosity, and this is taken into account in the model. But the experiments conducted earlier to confirm this point of 
view revealed a correlation between these indicators at the level of 6 % — and this is approximately equivalent to the 
percentage of errors made when using a mathematical model (in comparison with the calculated data). 

To express the cumulative effect of the fluid density and viscosity on the hydraulic brake operation, you can 
enter the parameter up in the model. However, this is not necessary, given the relatively low error rate. Otherwise, the 
proposed model will become much more complicated. 

Simulation of field conditions. The proposed mathematical model expresses a physical phenomenon 
identified and studied under the laboratory conditions at the initial braking torque 355 N-m. 

In the field, ballasts are characterized by an initial braking torque value Mg. It is proposed to introduce into 


the mathematical model a parameter that will not change its shape, but, presumably, will reflect the “field”: 

M, =58.5P + 284.69 +1.291+10.62+(M,, —M,, ) : (19) 
where Mz — the braking torque of the “field” ballast after applying the conditions (pressure, density, and viscosity); 
M po — the initial braking torque for the field ballast; /;,, — the initial braking torque of the laboratory ballast, equal to 
355 N-m. 


This is the main hypothesis. It is substantiated as follows. The dimensions of the laboratory inhibitor are not 
chosen randomly, but so as to correspond to the actual smaller dimensions of the brakes produced by Parmac L.L.C 


(model 112—500)'. 


' Original HYDROMATIC® Brake 620 251-5000. Parmac L.L.C // www.parmacbrake.com. 
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Discussion and Conclusions 

1. The mathematical model is created using the inverse analysis method, which relates the braking torque of 
the hydrodynamic brake to the operating parameters (pressure, viscosity, and density) of the fluid inside the ballast. 

2. The experiments with various random values of variables validated the adequacy of the model. The values of 
the braking torque determined experimentally and using the created model were compared. The error rate did not 
exceed 6 % (Table 6). 

3. An amendment to the model for simulating the operation of the hydrodynamic brake in the field 1s proposed. 

4. Based on the results obtained, it 1s advisable to test the created model in the field with a real payload in 
future studies. 
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Simulation of the hydraulic system of a device with self-adaptation 
for power and kinematic parameters on the working body 
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Introduction. Currently, Russia has adopted a course towards the creation of intelligent machines and equipment. The 
same holds for mobile technological machines for road construction and public utilities. Therefore, the design and 
creation of this type of actuators with a self-adaptation function is a critical task. 

Materials and Methods. A device equipped with a hydraulic drive with self-adaptation to load and coordination of 
kinematic and power parameters of the principal motion and the feed movement of the working body of the rock- 
drilling rig, is presented. To study and design the device based on the mathematical modeling methods of a hydraulic 
drive and adaptive systems, a mathematical model is proposed. It is developed using the foundations of the theory of 
volumetric stiffness of hydraulic systems. This enables to accurately describe the impact of the dynamic properties of 
the hydraulic system (compressibility of the working fluid, elastic properties of pipelines, high-pressure hoses, 
hydraulic apparatuses) on the dynamic properties of the system as a whole. 

Results. The mathematical model for a device with self-adaptation includes submodels of adaptive communication, 
interrelations of power, kinematic and process parameters of rock drilling, as well as mathematical description of the 
movement of system elements. The solution to the developed mathematical model was performed in the software 
environment for dynamic modeling of technical systems SimInTech. As a result, general dependences of the adaptive 
system on the design parameters of the system and the operating conditions are obtained. 

Discussion and Conclusion. The mathematical model of the presented device shows the fundamental possibility of 
implementing the principle of self-adaptation in terms of load under external and internal disturbing actions during 
operation. The results obtained can be used under designing adaptive systems of other technological equipment, for 
example, for the implementation of deep drilling in workpieces with variable properties in its depth. 


Keywords: hydraulic device with self-adaptation, hydraulic drive function, generalized mathematical model, adaptive 
communication, coordinated movements, working body, load, stabilization. 
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Introduction. The course adopted in Russia on building intelligent machines is the basis for the design and 
creation of actuators with a self-adapting function. The solution to this problem is timely and relevant. 

Drives with differential couplings of internal structure elements possess self-adapting properties [1]. This class 
of technical systems includes a variable rock drilling device (RF Patent No. 2582691). The self-adapting property is 
implemented by a device with negative feedback and positive feedback [2]. 

The quality of the self-adapting process is affected by external and internal actions — load variability, medium 
resistance, dry and viscous friction, volumetric stiffness of the fluid and pipelines, adaptive links. 

Materials and Methods. At the current level of development of computer technology, the complexity and 
high cost of the designed equipment, methods of mathematical modeling are widely used at the development stage [3- 
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13]. The development of a special model for calculations and computational experiments to determine the parameters at 
the design stage provides a reasonable choice of the standard size of the device with self-adaptation. 

The authors have developed a mathematical model for solving the problem of device design in the environment 
of dynamic simulation of SimInTech technical systems (Simulation In Technic) [14, 15]. This software product enables 
to simulate technological processes occurring in various industries with simultaneous modeling of control systems. 

Fig. 1 shows a device for exercising the function of self-adaptation under the conditions of force resistance on 
the working body which contains fixed displacement pump (H), security valve (KII), filter (M), adjustable throttles 
(J[p1) and (JIp2), flow regulator (PP), feed hydraulic cylinder (T1]) and main drive hydraulic motor (ITM), hydraulic 
valves (P1), (P2) and (P3), pressure gauges (Mu1), (Mx2) and (Mx3), 
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Fig. 1. Schematic diagram of the device for exercising the function of self-adaptation under the conditions of 
force resistance on the working body 
The device is designed to exercise the function of self-adaptive load (to stabilize it) and matching output 
movements working on production machines and equipment, as well as functions of the actuator itself. Therefore, the 
device provides multi-position distributors for the formation of flows and directions of the working fluid. Such 
additional device functions include heating of the hydraulic oil in its start operation; “weighting” of moving parts (idle 
load) of the feeder when adjusting to technological mode; accelerated lifting (retraction) of the tool with rotation, but 
without controlling its speeds; tool feed “to the object of action” or “from the object of action” with rotation and control 


of tool speeds. 
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The mathematical dependence for the internal negative coupling of the device is established through joint 
solving the equations arising from the pressure balance equations in the system and the working fluid flow balance 
equation in the feedback system. 


From the pressure balance in the system, it follows 





_ 1 1 frat 
AP aps —_ wn, em + Pa ~~ (1 > fi ) APypi0. (1) 
where Ap,»; and Ap,,19 — differential pressure across throttles 5 and 10; w,, and M,, — reduced volume of the 


hydraulic motor IM and the torque generated by it; f, and f,,,,, — areas of the hydraulic cylinder piston in the piston 
and rod cavities; F, — resistance force from the ground, overcome when moving the rod of the hydraulic cylinder 
ru. 

The balance equation for the costs in the feedback loop has the form 

Qnp1o = Qaps + Vn’ fare (2) 
where Q,p5 and Q,519 — the flow rate of the working fluid through the throttles J[p1 and J[p2; v, — the hydraulic 
cylinder I'L] piston travel speed. 

Considering that the hydraulic pump H, which feeds the system, has a constant capacity, the flow of the 
working fluid through the flow regulator PP under the operation of the drive system remains constant. It becomes 
obvious that the flow rate of the working fluid through the throttle J[p2 will also be constant. In this case, the 
dependence of the acceleration of the hydraulic cylinder piston I'L] on the total load on the hydraulic motors of the 


main movement and the feed movement will take the form 


dvy 1 dM , 1 dF 
‘dt — A (— dt - Sa Ea , ©) 





where v, — hydraulic cylinder II] piston travel speed; t — time; A — feedback coefficient depending on the 
parameters of the nominal operating mode of the system, the design parameters of hydraulic machines, throttles and 
its settings. 

It can be seen from the equation (3) that the positive increments of the torque on the hydraulic motor shaft 
and the movement of the hydraulic cylinder feed piston correspond to the negative increment of the tool feed speed. 
In other words, as the ground resistance increases to the rotation of the cutting tool or the movement of its feed, the 
tool feed rate decreases. 

The direct relationship between the speeds and accelerations of the hydraulic motor shaft of the main 


movement and the hydraulic cylinder rod of the feed movement is as follows 








1 
Y= oe (Wi Orn _ ‘cnn 
fh 
Or 
Avy AWry 
dt = Wry dt (4) 


where W,, — the angular speed of rotation of the hydraulic motor shaft; Q,, — the flow rate of the working fluid 
through the flow regulator. 

From the equation (4), it can be seen that with an increase in the rotation speed of the hydraulic motor of the 
main movement (occurs with a decrease in resistance from the side of the treated surface), the speed of movement of the 
feed cylinder piston increases. 

Mathematical modeling of the movements of the working elements of the device. The equations of motion 
of the rotor of the hydraulic motor of the main movement and the piston of the hydraulic cylinder of the tool feed have 


the form: 


AW ru 


Jom dt — APrm Wry — Meonp » (5) 
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where J,,, — the total moment of inertia of all rotating elements of the system reduced to the shaft of the hydraulic 
motor; m, — the mass of all moving parts of the system reduced to the hydraulic cylinder piston; M.onp) — the total 
moment of resistance to the rotation of the working body reduced to the shaft of the hydraulic motor; F.o,, — total 
force of resistance to the tool movement from the side of the treated surface reduced to the piston of the hydraulic 
cylinder; Ap,,,— pressure drop across the hydraulic motor; p, and p,,,,— pressure in the piston and in the rod 
cavities of the hydraulic cylinder, respectively. 

Modeling the properties of a hydraulic system 

The mathematical model of the hydraulic system of the device for performing the function of self-adaptation 
under the conditions of force resistance on the working body is developed using volumetric rigidity, which provides 
modeling as close as possible to the real characteristics [16—20]. When modeling, special attention is paid to 
determining the reduced coefficient of volumetric stiffness of high-pressure hoses [16]. The resulting mathematical 
model includes the following equations: 

— the equation of the pressure increment at various points (in Fig. 1, marked by points 1 to 23) of the 
hydraulic system has the form: 

p= Ci), Cat = 2s Cre 0; (7) 

where )) Q,,; and }) Qucxij —- the sums of all the flow rates of the working fluid entering and outgoing from the 
considered (i-th) volume of the system during the time dt; C,,; — the reduced coefficient of volumetric rigidity of the 
selected section of the hydraulic system; 


— the equation for determining the flow rate of the working fluid through various elements of the hydraulic 


2 : 
Q; = uf F Pi — Pi+il * Sign(p; — Di+1), (8) 


where p; and p;,; — pressure at the inlet and outlet of the hydraulic resistances; f — free cross-sectional area of 


system has the form: 


resistance; p — working fluid density; 


— the formula for calculating the reduced flow rate of linear resistances is as follows 


p= = ina (9) 
where d; and /,; — diameters and lengths of the linear section of the pipeline; A, — the coefficient of hydraulic friction of 
the pipeline; 

— the formula for determining the reduced coefficient of volumetric stiffness of the metal sections of the 


pipeline has the form: 








a ee 
C; _ td2 esl (10) 
5 E| 


where d and / — the internal diameter and length of the pipeline section; 6 — its wall thickness; Ey, and £;— moduli of 
elasticity of fluid and wall material of the hydraulic line. 
The reduced coefficients of the volume stiffness PBJ] are determined experimentally. 


The performance of the hydraulic pump 1s determined with account for the volume losses from the formula: 
dpuOu 
Qua a Now (11) 
where qg,, — the working volume of the hydraulic pump; ,, — the hydraulic pump shaft speed; no, — the current value 


of the volumetric efficiency of the hydraulic pump. 


The flow rate of the working fluid through the hydraulic motor is determined from the formula: 
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27 0,m 
where gp, — the hydraulic pump working volume; ,, — the hydraulic pump shaft speed; no, — current volumetric 
efficiency of the hydraulic motor. 
The current value of the volumetric efficiency coefficients of the hydraulic pump and the hydraulic motor are 


determined from the formula: 





No =1-(1—Nonom) *— > (13) 


Pp.nom 
where 10.nom — the nominal volumetric efficiency of the hydraulic pump and hydraulic motor; P,.nom — the value of the 
nominal pressure of hydraulic machines; p, — the current value of the pressure on the pump or motor. 

Research Results 

The proposed mathematical model of the device under consideration provides high-precision theoretical 
studies on the operational capabilities of rock drilling at the design stage. The calculation of the drilling system, 


performed using the SimInTech software [10, 11], showed the correctness of this statement. 
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Fig. 2. Changing the parameters of the drilling machine operation with a stepwise change in the resistance to rotation of the main 
working body of the drilling machine from the ground: a) — movement of the hydraulic cylinder piston (x,,) and its speed (v,); 
b) — torque of resistance to rotation of the hydraulic motor shaft M,,, 


Fig. 2 shows the graphs of change in the parameters of the piston movement of the feed hydraulic cylinder with 
a stepwise change in the moment of resistance to rotation of the main working body of the drilling machine from the 
ground. 

The analysis of calculation results presented in Fig. 2 shows that with an increase in the moment of resistance 
to rotation of the working body of the drilling machine (Fig. 2 5), the speed of translational movement of the tool 
decreases and can take negative values (Fig. 2 a). This indicates that the system performs the adaptation function. The 
properties of the system require additional research. 

Conclusion. The proposed mathematical model enables to make a preliminary performance assessment at the 
design stage and select the most rational parameters of the drilling system equipped with an adaptive hydro-mechanical 
drive under different operating conditions and varying its design properties. The results obtained and their analysis 
suggest that the proposed method of adapting the main and auxiliary movements of the drive system can be successfully 
applied in other technological equipment, for example, in deep drilling of multilayer metal workpieces. 

References 

1. Ivanov KS, Yaroslavtseva EK. Teoriya mekhanizmov s dvumya stepenyami svobody [Theory of 
mechanisms with two degrees of freedom]. Lambert Academic Publishing; 2014. 180 p. (In Russ.) 

2. Khinikadze TA. Issledovanie sootvetstviya kharakteristik adaptivnogo modulya_ gidroprivoda 
tekhnologicheskim parametram mashin [Study of characteristics conformity of the hydraulic drive adaptive module to 
technological parameters of the equipment]. Young Researcher of the Don. 2018;2(11):107—112. (In Russ.) 

3. Popikov PI, Goncharov PE, Sharov AV. Matematicheskaya model’ rabochego protsessa lesnogo 
pozharnogo gruntometa s ehnergosberegayushchim gidroprivodom [Modelling workflow forest fire soil-thrower with 
energy-saving hydraulic drive]. Forestry Engineering Journal. 2017;4:182—189. (In Russ.) 


Machine building and machine science 


7 
\O 


http://vestnik-donstu.ru 


Advanced Engineering Research 2021. V. 21, no. 1. P. 55—61. ISSN 2687-1653 





4. Bartenev IM, Popikov PI, Sharov AV, et al. Osobennosti chislennogo integrirovaniya sistemy 
differentsial'nykh uravneni rabochego protsessa pochvoobrabatyvayushchei mashiny s gidroprivodom dlya profilaktiki 
i tusheniya lesnykh pozharov [Features of numerical integration of the system of different equations of the working 
process of soil-cultivating machine with a hydraulic drive for the prevention and suppression of forest fires]. Forestry 
Engineering Journal. 2018;3(31):170—176. (In Russ.) 

5. Popikov PI, Scheblykin PN, Sharov AV, et al. Modelirovanie srabatyvaniya predokhranitel'nykh ustroistv 
lesnogo pozharnogo gruntometa [Modeling actuation safety devices forest fire soil - thrower]. Aktual'nye napravleniya 
nauchnykh issledovanii X XI veka: teoriya 1 praktika. 2017;5(1):336—340. (In Russ.) 

6. Sova AN, Stepanov MI, Sova VA, et al. Teoreticheskie osnovy modelirovaniya dinamicheskikh 
kharakteristik privodov perspektivnykh kosmicheskikh apparatov s uchetom funktsionirovaniya ostsillyatorov 
[Background for modeling the dynamic characteristics of advanced spacecraft drives considering the operation of 
oscillators]. Vestnik of DSTU. 2019;19(4):317-327. URL: https://www.elibrary.ru/item.asp?id=41591369 (accessed: 
27.01.2021). (In Russ.) 

7. Syrkin VV, Kvasov IN, Galuza YuF, et al. Issledovanie dinamiki delitelya potoka s reguliruyushchim 
organom iz ehlastomera [Study of influence of fluid compressibility on dynamics of flow divider]. Omsk Scientific 
Bulletin. Series “Aviation-Rocket and Power Engineering”. 2018;22(4):9-14. (In Russ.) 

8. Alekseev GV, Aksenova OI. Vozmozhnosti modelirovaniya oborudovanitya dlya_ snizheniya 
ehnergoemkosti realizuemykh tekhnologicheskikh protsessov [Equipment modeling to reduce technological processes 
power consumption]. Journal of International Academy of Refrigeration. 2015;4:49—54. (In Russ.) 

9. Berdnik Y, Beskopylny A. The approximation method in the problem on a flow of viscous fluid around a 
thin plant. Aircraft Engineering and Aerospace Technology. 2019;91(6):807—813. 

10. Beskopylny A, Kadomtseva E, Strelnikov G, et al. Model of heterogeneous reinforced fiber form concrete 
in bending. In: IOP Conference Series on Materials Science and Engineering. 2018;365(3):032023. 

11. Beskopylny AN, Kadomtseva EE, Strelnikov GP. The boundary condition influence on a stress-strain state 
of a corrugated plate on an elastic foundation. Materials Science Forum. 2018;931:60—65. 

12. Beskopylnyi AN, Lyapin AA. About the identification of layered constructions properties. Recent Trends 
in Science and Technology Management. 2016;2:36—45. 

13. Chukarin AN, Beskopylny AN, Isaev AG. Issledovaniya vibroakusticheskikh kharakteristik v rabochei 
zone operatora pri abrazivnoi obrabotke svarnykh shvov [Research of vibroacoustic characteristics in the operator 
working area during welds abrasive treatment]. Occupational Safety in Industry. 2019;11:7—12. (In Russ.) 

14. Kartashov BA, Kozlov OS, Shabaev EA. Dynamic mode environment of technical systems SimInTech: 
manual. Moscow: DMK-Press; 2017. 424 p. 

15. Abalov AA, Nosachev SV, Zharov VP, et al. Using the SimInTech dynamic modeling environment to build 
and check the operation of automation systems. MATEC Web of Conferences. 2018;226:04003. 

16. Rybak AT, Shishkarev MP, Demyanov AA, et al. Modeling and calculation of hydromechanical systems 
dynamics based on the volume rigidity theory. MATEC Web of Conferences. 2018;226:01001. 

17.Rybak AT, Temirkanov AR, Lyakhnitskaya OV. Dynamics of Synchronous Hydromechanical Drive in 
Mobile Machine. Russian Engineering Research. 2018;38(9):702—704. DOI 10.3103/S1068798X 18090253 

18. Temirkanov AR, Rybak AT. Sinkhronnyi gidromekhanicheskii privod rabochego organa mobil'noi mashiny 
1 ego matematicheskaya model' [Synchronous hydromechanical drive of mobile actuator and its mathematical model]. 
Vestnik of DSTU. 2013;13(3-4):104—-110. URL: https://www.vestnik-donstu.ru/jour/article/view/400 (accessed: 
27.01.2021). (In Russ.) 

19. Temirkanov AR, Rybak AT. Modelirovanie sistemy gidromekhanicheskogo privoda rabochego organa 
mobil'noi tekhnologicheskoi mashiny [Modeling of hydromechanical drive system of working body for mobule 
technology machine]. Vestnik of DSTU. 2014;14(4):176-185. URL: https://www.vestnik- 
donstu.ru/jour/article/view/362 (accessed: 27.01.2021). (In Russ.) 

20. Mirny VI. Nekotorye voprosy povysheniya proizvoditel'nosti gidravlicheskogo privoda perforatsionnogo 
pressa tsiklicheskogo deistviya [Several issues of productivity increase and economic effectiveness of hydraulic 
perforated hammer of cycling capacity]. Vestnik of DSTU. 2008;8(1):75-85. URL: https://www.vestnik- 
donstu.ru/jour/article/view/1212 (accessed: 27.01.2021). (In Russ.) 


Submitted 04.12.2020 
Scheduled in the issue 21.01.2021 


Khinikadze T. A., et al. Simulation of the hydraulic system of a device with self-adaptation for power and kinematic parameters 





About the Authors: 
Khinikadze, Tengiz A., postgraduate student of the Instrument Making and Biomedical Engineering 
Department, Don State Technical University (1, Gagarin sq., Rostov-on-Don, 344003, RF), ORCID: 


https://orcid.org/0000-0003-1709-9505, khinikadze@mail.ru 


Rybak, Alexandr T., professor of the Instrument Making and Biomedical Engineering Department, Head of 
the Certification and Professional Accreditation of Educational Programs Department, Don State Technical University 
(1, Gagarin sq., Rostov-on-Don, 344003, RF), Dr.Sci. (Eng.), professor, ORCID: https://orcid.org/0000-0001-9950- 


3377, 2130373(@mail.ru 


Popikov, Petr I., professor of the Forestry Mechanization and Machine Design Department, Voronezh State 
University of Forestry and Technologies named after G.F. Morozov (8, Timiryazeva St., Voronezh, 394087, RF), 
Dr.Sci. (Eng.), professor, ResearcherID: AAX-7892-2020, ORCID: _ https://orcid.org/0000-0002-6348-8934, 


popikovpetr@yandex.ru 


Claimed contributorship 

T. A. Khinikadze: basic concept formulation; research objectives and tasks; computational analysis; text 
preparation; formulation of conclusions. A. T. Rybak: academic advising; analysis of the research results; the text 
revision; correction of the conclusions. P. I. Popikov: analysis of literature and patent materials, dynamic characteristics 
of the hydraulic drive. 


All authors have read and approved the final manuscript. 


Machine building and machine science 


ON 
— 


http://vestnik-donstu.ru 


Advanced Engineering Research 2021. V. 21, no. 1. P. 62—70. ISSN 2687-1653 





UDC 621.791.75:004.942 https://do1.org/10.23947/2687-1653-202 1-2 1-1-62-70 


Consumable electrode — additional filler wire arc interaction control under 
surfacing (DE-GMAW) 


S.S. Poloskov’, V.A. Erofeev’, M. A. Sholokhov" 


‘Moscow Institute of Physics and Technology (National Research University) (Moscow, Russian Federation) 





* Tula State University (Tula, Russian Federation) 


> Ural Federal University (Ekaterinburg, Russian Federation) 


Introduction. Arc surfacing through feeding an additional filler wire heated by an additional arc burning between the 
filler wire and the electrode wire is considered. Under the conditions of such surfacing, the minimization of the input of 
the remelted substrate metal into the weld metal is studied. The research objectives are to examine the conditions 
providing self-regulation of this advanced arc surfacing process, and to evaluate control capabilities of the heat impact 
power on the metal and on the weld metal flow. 

Materials and Methods. In solving a wide range of welding and surfacing tasks, it is advisable to use engineering 
analysis methods based on physicomathematical modeling of processes and phenomena. These include: 

— self-regulation of the arc process under joint melting of the electrode and the filler wires; 

— assessment of the possibilities to control the heat impact power on the metal and on the weld metal flow during the 
formation of the weld pool. The features of the arc surfacing of anticorrosive chromium-nickel steels on low-alloy steel 
are considered in the paper. 

Results. New mathematical dependences are proposed that describe physical phenomena under surfacing with an arc 
interaction between the electrode and filler wire. A physicomathematical model of the joint melting of the electrode and 
filler wire is developed. It provides determining the values of the control parameters. In addition, you can find out how 
much heat affects the substrate from: 

— heat release in the main arc, 

— droplet flows of the weld electrode and filler metal, 

— arc plasma radiation. 

Discussion and Conclusions. It 1s established how the current and the lengths of the main and additional arcs are 
affected by the supply voltages. The feed rate of the electrode and filler wire with a diameter of 1.6 mm and 1.2 mm 
made of Inconel 625 alloy is determined. It is shown what thermal effect the substrate undergoes in this case. It is noted 
that due to the larger value of the main arc current, the diameter of the electrode wire should be larger than that of the 


filler wire. The heat flow in the substrate is created mainly by the flow of the weld metal droplets. 
Keywords: physicomathematical model, arc interaction, surfacing, consumable electrode, filler wire. 
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Introduction. Surfacing with materials whose characteristics differ markedly from the base metal, provides the 
necessary performance properties of the contact surfaces of the products. When surfacing, it is required to reduce the 
mixing of materials, so you need to adjust the power of the heat impact on the surface of the base metal, and the amount 
of melting material. In [1], the possibility of such regulation through changing the heat flow during the arc interaction of 
a melting electrode with an additional filler wire was first validated. Initially, this approach was used mainly in welding 
[2-4], and in recent years — in surfacing [5]. According to the author of the paper [6], this method of surfacing will help 
to provide high performance properties of the sealing and contact surfaces of pipeline fittings. However, a number of 
researchers have noted a factor that complicates the use of this method in welding and surfacing. They consider the 
major problem to be the filler wire feed control. At an insufficient feed rate, it can melt outside the molten pool, and at a 
high feed rate, it cannot melt. Therefore, to provide the required flow of the deposited material with a limited heat 
impact of the arc on the metal surface, it is necessary to accurately set the current of the additional arc burning between 
the wires. The working conditions, the instability of the wire feed, and the heterogeneity of their thermophysical 
properties also unpredictably change the melting conditions of the main electrode and the filler wire. Therefore, the 
surfacing process stability can be provided through using the self-regulation effect of the main and auxiliary arcs. The 
work objectives are to study conditions for the self-regulation of the arc process through co-melting of the electrode and 
filler wire, as well as to assess the possibility of heat power control for the metal and the filler metal flow. 

Materials and Methods. When solving a number of problems, the use of engineering analysis methods based 
on physical and mathematical modeling of welding processes is challenging [7, 8]. Such models are a system of 
differential equations, whose boundary conditions take into account many technological factors. The equations are 
solved in an iterative cycle. An important advantage of this method of engineering analysis based on fundamental 
physical laws is the universality of the results and the possibility of using them to study the mechanisms of physical 
interactions in welding and surfacing processes [9]. Such phenomena include: 

— self-regulation of the arc process under the joint melting of the electrode and filler wires, 

— assessment of the possibility of controlling the power of thermal impact on the metal. 

However, taking into account the effect of arc self-regulation, it is required to further investigate features of 
the joint melting of the electrode and filler wires. In addition, you need to know with what power the substrate is 
affected by heat from: 

— heat generation in the main arc, 

— droplet flows of the deposited electrode and filler metal, 

— arc plasma radiation. 

In this paper, we consider double-electrode gas metal arc welding, DE-GMAW, and electrode surface welding 
(a melting electrode and a current-carrying filler wire) in a protective gas. 

Main Part. To study the mechanism of interaction between the main and additional arc between the electrode 
and filler wires, a physicomathematical model is required with account for the following: 

— features of the melting of the electrode and filler wires, 

— the effect of self-regulation of the thermal power of arcs under surfacing. 

The numerical solution to the model equations in an iterative cycle provides a deeper study of the adjustment 
characteristics of the process and their influence on the thermal power of the surfacing. 

Physical model of electrode and filler wire melting. Taking into account many significant factors, it is 
appropriate to use the approaches described in [10] when studying the features of self-regulation of melting of electrode 
and filler wires. However, the modeling space needs to be specified. Consider the refined modeling space for the 
melting conditions of the electrode and filler wire under the action of arcs. The first of them — a (Fig. 1) burns between 
the electrode wire | and the substrate 3, and the second b — between the electrode and the filler wire 2. In this case, in 
our opinion, it is suitable to feed the wire from the front of the molten pool. This will eliminate welding (“freezing’’) of 


the unfused filler wire in the solidifying metal of the tail section of the molten pool. 
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Fig. 1. Process chart of double-arc consumable electrode surfacing with an additional filler wire 


To form a common plasma column of arcs a and 5, the filler wire is fed at an angle o to the axis of the melting 
electrode x. This enables the axes of the filler and electrode wires to intersect at a given distance X from their current 
contact nozzles, and to set the distance x, from the intersection point of the wire axes to the metal surface. In the 
simulation space, feed rates and wire diameters of the electrode wire v,, d, and the filler wire v, dy are also set. 

For the main arc a, voltage U-, is given, and for the additional arc b — U,,. These voltages are regulated by the 
key K, which redistributes voltage U, of the power supply between the arcs. The continuity of arc currents J, J; are 
provided by chokes L,, Ly and closing diodes D,., D; Arc currents /,, Jy are self-adapted through a self-regulatory effect 
of lengths /., J; of the corresponding arcs. The lengths of the electrode stickouts X,, X; are determined through the 
melting of the electrodes. The melting rates are determined by the arc currents /,, J; and the voltage drop in the anode 
regions, as well as by the heating of the stickouts by the arc currents passing through them. R — electrical resistance of 
the power supply. 

Mathematical model of the surfacing process self-regulation. The main control action on the process is 
provided by the operation of the key K. Its closures (with a period of t and duration of ¢,, ¢) regulate the voltage on the 
electrodes: 

Usp = (Ue — (Ret R)le) £5 User = (Ue — (Rp + Ret RY) L, (1) 
where R,, Ry— electrical resistances of the electrode and filler wire stickouts. 

In the steady-state mode, these voltages determine the arc lengths [11]: 

le = = (Usp — Ua — Ux): lp = = (Uer — Ua — Ur), (2) 
where U,, U;, — voltage drop in the anode and cathode regions of the arc; E — potential gradient in the arc column. 

It is advisable to determine the average values of the arc currents through the melting rates of the wires, which, 
in the steady-state mode, are equal to the rates of their supply: 


2 
Vemde 


nd? 
Ip = (co(Te — Tay) + Ays)ste = “Eh * (OTe — Tre) + 415); (3) 








U 
where 7;, — temperature of electrode metal drops; cp — heat capacity of metal of wires; g,s — melting heat; Ty, Ty; — 
preheat temperatures of wire stickouts by arc currents. 

It should be taken into account that with pulse regulation, the current will change with the period t and the 
amplitude: 


A= =. (4) 
However, with a small period of change, typical for modern power supplies, the periodic change in current can 
be ignored. 
Current flowing through the cathode spot on the metal surface: 
Le ale Ip. (5) 


To find out the temperature distribution along the stickout axes x of the electrode 7, (x) and filler 7;(x) wires, it 
is required to consider the dependence of the electric resistivity of the wire metal on the temperature p, (7) [12]. 
Then the temperature distribution 1s: 
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rae = (48) 


e 


(5 be (Tox) ax, Ty (2) = (“) ff pe(TeCe)) ae 6) 


The preheat temperature 77, is determined from its values at the end of the stickout x =X, 
As noted in [13, 14], the temperature of the metal drops 7; entering the weld (molten) pool is determined by a 


Vecp vecp © 0 


combination of many factors. However, in any case, it will be in the range between the melting point 7; and the boiling 
point 7) of the metal: 


Te (7) 
The heat dissipation power is calculated using the following formulas. 
— For the cathode arc spot: 
Px = Ure — Ip) (8) 
— Transferred by electrode droplets: 
nd? 
F, = Ugle+ cpl Ve he (9) 
— From the flow of filler wire droplets: 
ye 
The arc column radiation power: 
Prag = E(Iply + Tele). (11) 


Only a part of the arc radiation reaches the surface of the substrate. Note that the radiation intensity is inversely 
proportional to the square of the distance from the radiation center. Assume that this center is removed from the 
substrate surface by the length of the main arc. Then, the Powe on the surface of the substrate: 


Prad 
a st oe (12) 
where S — surface area, r — the distance from the sce to the arc torch axis. 








The numerical solution to the model equations is reduced to an iterative selection of arc currents, in which 
equality is achieved in all relations (Fig. 2). 


Source data: metal grade and surfaced material, 
dimensions X), X>, x. the surfacing mode parameters v,, d., va dy Ue, R, T, te, ty. 
Material properties from the databases 7), 77, cp, drs, Ua, Ux, E, pe (T). 
Initial approximation 
Iterative cycle of current refinement J,, J; 
Calculation of temperature distribution in stickouts: x = 0, T(x) = To, Re = 0, Ry= 0 
Cycle x =x + dx 
pe(T (x) 
vFCP 
Electric resistivity of the aaa 
De pe (7; (x)) 


Pe(Te(x) 
ee ls, sKA(x)= R(x) + ade 


dx 


T(x) =T a(x) + (Hs) pe) a T(x) = Ty) + (“4) 


RzQ)—hoG)) + 


as long as x < X, or x < X; 
Parameters of stickouts 7. = T.(L.), Tip = T(L)), Re = Ri (Le), Re = Ry (Lp) 
Clarification of currents 


1 VT Votd 
Ir —~ | 3Tp + AU, F (colt, — Ty) *4s)) le = 7 [ste +" + — 8 (ep Tx) +45) 


until /, const or J; const 
Determination of process parameters: D = kp - a i,)3 ras 
P, =Uale + cpTy,Ve EP = Uyly + ep yf sy 
P. =F (Iply + Tele) Wap be = | (Uep — Ya — Ur )i ly = 5 (Ver Ua ~ Ur) 


Fig. 2. Algorithm for calculating the parameters of the joint melting process of electrode wires 





The algorithm for calculating the parameters of the joint melting process of electrode wires presented in Fig. 2 
enables to study in detail the conditions of self-regulation of the main and additional arcs under surfacing. 
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Research Results 
Process control characteristics. Regulatory actions: 


— the feed rate of the electrode and filler wires, 
— the no-load voltage of the power supply, 


— the relative duration of the key closures, which controls the ratio of the voltages of the main and additional 
arcs. 

These actions determine the arc currents, their lengths, and the power distribution of the process. Additional 
parameters that determine the results of regulation are the wire diameters, their stickouts, and the angles between them. 

The ranges of control parameters values are limited. The most important of them is the condition for the 
existence of an arc between the substrate and the electrode, i.e., the current J, = J, —J;> 0 through the cathode spot on 
the substrate. 

Accordingly, the relation should be fulfilled: 

ved. Urey (13) 
Ua Uk 

The second critical limitation is the minimum voltage on the arcs, which provides their stable arcing in the 

absence of short circuits: 


Ue “2 > Ug + Ups Ue £ > Unt Up. (14) 
In addition, it 1s required to provide the correct location of the vanishing point of the electrodes above the 
substrate surface and the deposited layer, and the absence of short circuits by the electrode metal drops, 1.e., the length 
of the main arc should be: 
1, >xXe+ ipctg—; Xp >do. (15) 
To do this, the main arc voltage should be at least: 
Ue © > Ugt Uyt Ele. (16) 


The length of the additional arc /; is a process parameter that provides removing the flow of additive drops 
from the main arc. 

Influence of the control parameters on the process power. The power of the heat fluxes acting on the substrate 
in the sum of the heat release in the cathode spot of the main arc, the radiation power of the main and additional arcs, 
heat transfer by drops of the electrode and filler metal. Fig. 3 shows the effect of the feed rate of the electrode wire with 
a diameter of 1.6 mm made of Inconel 625 alloy on: 

— current /, of the main arc, 

— power P,, emitted in the cathode spot, 

— power P, of the heat flux drops of electrode metal, 

— temperature T, of heating the electrode stickout, 

— length /, of the main arc. 
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Fig. 3. Dependence of the main arc parameters on the feed rate v, of the electrode with a diameter of d, = 1.6 mm made of 
Inconel 625 alloy at stickout of /, = 10 mm, supply voltage U, = 32 W, duty cycle ¢,/t = 0.5, and the additional arc current J-= 108 A 
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With an increase in the speed v, of the electrode wire feed at a fixed current /; of the additional arc, the current 
I., flowing through the electrode, increases almost in proportion to the speed. The deviation from the proportionality is 
explained by a noticeable increase in the temperature T, of the electrode stickout. For the same reason, the power P, that 
heats the electrode increases faster than the arc current J,. The power P;, released in the cathode spot on the substrate 
surface is almost proportional to the feed rate, but it is much less than the thermal power P, of the droplet flow. The 
length /, of the main arc decreases linearly as the electrode feed rate increases. 

Fig. 4 shows the impact of the feed speed v, of the filler wire with a diameter of 1.2 mm made of Inconel 625 
alloy on: 

— additional arc current J;, 

— power P, generated in the cathode spot, 

— power P; of the heat flux drops of metal filler wire, 

— temperature 7; of the heated wire stickout, 


— length /, of the additional arc. 
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Fig. 4. Dependence of the additional arc parameters on the feed speed v, of the filler wire with a diameter of d, = 1.2 mm made of 
Inconel 625 alloy at stickout of /, = 10 mm, supply voltage U, = 32 W, duty cycle t/t = 0.5, and the main 


With an increase in the feed rate v, of the filler wire at a fixed current /, of the main arc, the current /; is 
somewhat disproportionate to the feed rate due to a noticeable increase in the temperature 7; of the wire stickout. The 
power P., spent on heating the wire increases faster than the current /; of the additional arc. The power P;, generated in 
the cathode spot on the substrate surface, is greatly reduced through increasing the feed rate of the filler wire due to the 
redistribution of the main arc current to the filler wire. The thermal power P; of the flow of the droplets from the filler 
wire is less than the power P, of droplets from the electrode. The length /; of the additional arc decreases non-linearly as 
the feed rate increases. This is due to the impact of the electrical resistance, which is greater in the filler wire than in the 
electrode stickout due to the smaller wire diameter and higher temperature 7; The voltage of the arcs determines their 
length, 1.e., their location above the substrate. 

Fig. 5 shows how the distance between the electrode and filler wire and their distance from the substrate 


depend on the surfacing mode. 
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Fig. 5. Impact of the main arc voltage U,, at different additional arc voltage U., on: arc radiation power P,qq , arc radiation power on 
the substrate P,,,, main arc length /,, main arc length x, from the vanishing point of the electrodes to the substrate 

The minimum arc voltage (12 V) is determined from the sum of the anode and cathode voltages, as well as 
their minimum length. The voltage U,, is set between the electrode wire and the filler wire. At the same time, an 
increase in the voltage U., between the electrode and the substrate increases the length /, of the main arc, which is 
accompanied by an increase in the power P,q of heat release in the plasma arc torch. Rising the voltage U., between the 
electrode and the additive increases the length /; of the additional arc, which also increases the power P,ag. Some part P,, 
of this power P,,q 1S radiated to the substrate surface. The dependence of the total power P,,q of the arc torch on the 
voltage of the main arc U,, is noticeable. At the same time, a change in the power P,,, that heats the substrate is less 
significant. This is due to the removal of the arc plasma from the substrate as the length /, of the main arc increases. The 
voltage of the additional arc U., affects the heating power P,,, of the radiation much more strongly, since the distance 
between the arc torch and the substrate does not depend on this voltage. Voltages U., and U-., affect significantly the 
location of the wire vanishing point x,, which imposes restrictions on the power supply of the arcs. 

Discussion and Conclusions. Surfacing with the electrode and filler wire arc interaction enables to adjust the 
power consumed for melting the deposited metal and for heating the substrate. The stability of consumable-electrode 
arc welding is usually provided by the self-regulating effect, which equalizes the melting and feeding rates of the 
electrode wires when they change over a wide range. The physicomathematical model of joint melting of electrode and 
filler wires is developed with account for self-regulation, in which the control actions are the wire feed rates and the 
voltage on the welding torches, while the currents and arc lengths are the results of self-regulation. The model provides 
determining the values of the control parameters, as well as the power of the thermal effect on the substrate under the 
consumable-electrode surfacing with additional filler wire. 

Based on the results of the work performed, the following conclusions can be drawn. 

1. Under the consumable-electrode arc surfacing with additional filler wire, to effectively regulate the ratio 
between the volumes of the surfaced material and the remelted substrate, the following self-regulation should be 
provided: 

— of the main arc between the electrode and the substrate, 

— of the additional arc between the electrodes. 

2. Self-regulation of the arc melting process of the electrode and filler wires is described by a system of 
equations. In it, the control actions are the wire feed rates and the voltages of the main and additional arcs, and the 
control results are the arc currents, their lengths, the distribution of heat generation power, the heat fluxes of the 


electrode metal droplets and the arc radiation on the substrate surface. 
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3. Limits are set on: 

— the ratio of the parameters of the feed of the electrode and filler wire, which provide the presence of an arc 
between the electrode and the substrate; 

— minimum voltage across arcs maintaining their stable burning in the absence of short circuits. 

4.Calculations of the modes of consumable-electrode surfacing with an additional filler wire of the Inconel 625 
alloy are performed. The calculations show that the process can be controlled over a wide range of wire feed rates and 


burner voltages. 
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Introduction. The relevance of the presented paper is due to the widespread use of plunger pumps in industrial practice, 
in particular, in gas and oil production. The quality of working operations and the efficiency of further well operation 
depend largely on their reliability. The improvement of plunger pumps involves increasing their reliability, increasing 
their service life, efficiency, downsizing, reduction in weight, labor intensity of installation and repair work. The 
modernization of the mechanism includes its power study since the found forces are used for subsequent strength 
calculations. Before the appearance of programs for the numerical analysis of solid objects, the analytical solution to the 
problem of strength calculation of the high-pressure pump drive frame was a very time-consuming and expensive 
procedure. The situation has changed with the development of computer technologies and the inclusion of the finite 
element method in the computer-aided design systems. The objective of this work is to perform a strength calculation 
on the TWS 600 plunger pump body made of 09G2S steel. 

Materials and Methods. A method for determining the reactions of the crank shaft supports of a high-pressure plunger 
pump and strength calculation of the drive part housing is developed. The direction and magnitude of the resulting 
forces and reactions of the supports are determined graphically according to the superposition principle of the force 
action on the supports. Strength calculations were performed using the finite element method in the computer-aided 
design system Solid Works Simulation. In this case, solid and finite-element models of the body with imposed boundary 
conditions were used, which were identified during the analysis of the design and the calculation of the forces arising 
under the pump operation. 

Results. The reactions in the crankshaft supports are described with account for the forces generated by the plunger 
depending on its operating mode and the crank position. The forces acting on each of the plungers and the resulting 
reactions in each of the supports are determined. The diagrams of stresses and the safety factor are presented, which 
provide assessing the strength of the body and developing recommendations for creating a more rational design. 
Discussion and Conclusions. As a result of the calculations, we have identified areas of the structure with minimum 
safety factors, and areas that are several times higher than the recommended values. This provides optimizing the design 
under study through strengthening the first and reducing the thickness of the metal on the second. From the point of 
view of weight and size characteristics and maintainability, the results of the strength calculation performed can be used 


to optimize the design of the pump body under typical operating conditions. 


Keywords: plunger pump, support reactions, strength calculation, design optimization, superposition principle, 
calculation of body parts. 
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Introduction. Oil production is not complete without the use of plunger pumps. They are required for such 
working operations as well cementing and acidizing, sand blast perforation, hydrofracturing, etc. [1]. Downsizing and 
reduction in weight of the pumps make them extremely attractive for use in mobile oilfield installations. The reliability 
of the pumps largely determines the quality of processing and the efficiency of further operation of oil and gas wells [2]. 

Despite the overall rather high level of plunger pump designs, they continue to be improved. The improvement 
of plunger pumps involves increasing their reliability, increasing their service life, efficiency, downsizing, reduction in 
weight, labor intensity of installation and repair work. For this purpose, the designs of the components and parts of the 
drive and hydraulic systems are changed [3-5]. The modernization of the mechanism includes its power study since the 
found forces are used for subsequent strength calculations. Recently, in connection with the development of computer 
technologies, numerical methods of strength analysis with the use of applied programs are increasingly used [6, 7]. 

Before the appearance of programs for the numerical analysis of solid objects, the analytical solution to the 
problem of strength calculation of the high-pressure pump drive frame was a very time-consuming and expensive 
procedure [8-10]. The development of the finite element method (FEM) in the deformable solid mechanics and its 
inclusion in computer-aided design systems (CAD, e.g., Solid Works Simulation) opens up new opportunities in solving 
problems of this kind. Strength calculations performed in CAD provide optimizing the design of the body parts. For the 
design under consideration, the body is one of the critical parts that takes up the load and provides the correct mutual 
arrangement of the drive part elements. 

To obtain reliable results in the strength analysis of the pump body using numerical methods in CAD, it is 
required to determine all external loads acting on the body. 

The widely used TWS 600 triplex plunger pump was selected as the object of modernization. In Russia, 
structural low-alloy 09G2S steel is used for manufacturing body parts of high-pressure plunger pumps [11]. Replacing 
the material of the body parts of TWS 600 plunger pump with 09G2S steel will reduce the cost. In this regard, this work 
objective is to analyze the possibility of replacing the material of the body of TWS 600 plunger pump with 09G2S steel 
to optimize the price and provide the possibility of repair work. 

Materials and Methods. The replacement of the material should be validated by the strength calculation. The 
calculation of the body was carried out by the FEM in the CAD Solid Works Simulation. 

The external forces acting on the pump body are reactions in the crankshaft supports. They arise from the 
action of the inertia forces of the reciprocating moving parts of the crosshead connecting rod group and the forces of 
fluid pressure on the plunger. To determine the reactions in the crankshaft supports, it is required to perform a dynamic 
calculation of the slide-crank mechanism of the plunger pump’. 

As mentioned in [12, 13], tasks of the dynamic analysis include studying the influence of external and internal 
forces on the links and kinematic pairs of the mechanism, as well as identifying ways to reduce dynamic loads. 

The crankshaft in question consists of two main journals on the major shaft axis and three crank journals 
eccentrically arranged at 120° offset. The bearing supports are located on the main journals and on the crank webs 


between the crank journals (Fig. 1). 


* Timofeev GA. Theory of machines and mechanisms: a course of lectures. Moscow, 2010. 351 p. (In Russ.) 
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Fig. 1. Crankshaft of TWS 600 plunger pump: | — main journals; 2 — crank journals; 3 — webs 


During the operation of the pump, three plungers reciprocate successively to provide discharge or suction (Fig. 2). 





Discharge 
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Fig. 2. Design diagram of the mechanism of the drive part of the plunger pump: a) position of the Ist plunger; 5) position of the 2nd 
plunger; c) position of the 3rd plunger 
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Consider the position of the mechanism in which the first plunger is located at the extreme discharge point (ex. 
p.) and is under maximum load. The angle of rotation of the crank journal for this plunger is zero (Fig. 2 a). Then, the 
second plunger performs suction (Fig. 2 5), its journal is displaced by 120°. The third plunger operates for the discharge 
(Fig. 2 c), the journal is displaced by 240°. Thus, in one situation, the crankshaft is loaded with forces of different 
magnitude, and the reactions of the supports will be different. This occurs when the crank is in the position when the 
first plunger is at the extreme point of discharge, the second performs suction, and the third operates for discharge. We 
will perform the calculation for the specified position of the crank (mechanism) [8-10], and then, for other positions: 
when the second and third plungers will reach the extreme discharge point. 

Link masses were calculated in the Autodesk Inventor Professional 2018 program based on the previously 
created 3D model of the drive of the TWS 600 plunger pump. A design diagram of the mechanism of the drive part of 
the plunger pump was constructed with an indication of the forces acting on its links (Fig. 2). 

Let us denote the external forces for finding the crank shaft reactions: P,, — fluid pressure force, P,, — inertia 
force. 

The fluid pressure force acting on the plunger is found from the expression: 

P= Pert * t, 
where p,,. — pressure of the pumped liquid applied to the plunger; r — radius of the section of the plunger’. 

The inertia force of the reciprocating moving parts 1s: 
P,= —m-i= —m-r-w*-(cos@ + (2: cos?@ — 1), 
m=m, + 0,275 -m,, 
where m; — crosshead mass, m,, — mass of the upper connecting rod head, r — crank radius, w — crankshaft speed at 
230 rpm, » — crankshaft angle. 
The force acting on the plunger is determined from the formula: 
Ps = Pai + Pui, 

where P,,; — pressure force of the pumped liquid in each position of the plunger, P,,; — inertia force of the 
reciprocating moving parts in each position of the plunger. 

Total force Ps, applied to the crosshead pin axis and directed along the axis of the cylinder, can be decomposed 
into: 

— force N acting perpendicular to the axis of the cylinder; 
— force S' acting on the axis of the connecting rod. 

N presses the crosshead against the cylinder wall, which causes wear on their surfaces. This force changes in 
value and direction, alternately pressing the crosshead to one or the other side of the guides. 

The force S, transferred to the of the crankpin axis, can be decomposed into: 

— tangential force 7, acting perpendicular to the crank of the crankshaft; 
— radial force K, directed along the crank axis. 

The force K 1s determined from the formula: 
cos(; + B;) 

Kj = Py; aT 
where ~; — crank angle (Fig. 2); 8; — angle of deviation of the connecting rod from the axis (Fig. 2); Ps; — total force 
acting on the plunger. 


The results of calculating radial forces for all crank journals in each considered position of the mechanism are 


summarized in Table 1. 


'Catalog of high-pressure plunger pumps manufactured by Weir SPM, Nord-SPM LLC. URL: http://suepcnm.pd/catalogues/JI2IK_ KATAJIOL 
HACOCOB.pdf (accessed: 14.02.2021). (In Russ.) 
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Table 1 
The results of calculating radial forces for all crank journals 


in each considered position of the mechanism 


~, degree 


oo [wo 
K,H 


The reactions of the supports are directed opposite to the radial forces (Fig. 3), and their values are determined 





from the formulas: 


_ 4 
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where K,, Kz, K3 — radial forces directed along the crank axis. 
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Fig. 3. Scheme of reactions of supports 


The direction and magnitude of the resulting forces F, and Fz of the reactions of the supports O,03 re 
determined graphically (see Fig. 3) according to the principle of superposition of the force action on the supports. 

The resulting forces F, and F, coincide in magnitude and direction with the forces R, and R3, respectively. 

Numerical values and directions of reference reactions are found for the position of the mechanism at which 
the first plunger was at the discharge extreme point (see Fig. 2 a), the second performed suction (see Fig. 2 5), and the 
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third — discharge (Fig. 2 5). For other positions of the mechanism, the reactions of the supports are determined 


similarly and will be numerically equal, alternately changing places. 
Table 2 shows the results of calculating the reactions of the crankshaft supports for the positions of the 


mechanism in which each of the plungers is alternately located at the extreme discharge point. 
Table 2 


Numerical values of the reactions of the supports for all positions of the mechanism 


Reaction forces of the supports 


Me anen Plunger position 
(2 
I ee : 225000 225000 153343 153343 
an 


II 


i. 2 
oe oe 153343 | 329602 | 225000 | 1251 
discharge suction 


1251 153343 329602 225000 


The body of the high-pressure plunger pump was modeled in Autodesk Inventor Professional 2018 CAD 
system. The model includes bearing cages since the presence of a cage significantly affected the nature of the load 





application. 
The strength analysis is performed by the FEM in the computer-aided design system Solid Works Simulation. 


Fig. 4 a shows a solid model of the calculated structure, and Fig. 4 b shows a finite-element model. 


Bearing cages 


Pump body 





b) 


Fig. 4. Body models of the TWS 600 three-plunger high-pressure pump: a) solid-state, b) finite-element 
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The finite element model contains 100,747 elements (175,172 nodes). The body material is 09G2S steel. The 
model describes the loads and records the movements in all directions of the surfaces of the bolt holes. Fig. 5 shows the 
boundary conditions for the contacting surfaces. By the type of contacts, the surfaces are connected to common nodes. 








6 


Fig. 5. Boundary conditions: 1 — operating load (radial force distributed over cylindrical surfaces); 2 — remote loading 
(mass of the gearbox and hydraulic unit); 3 — restriction — “fixed”, zero movements in all axes —X, Y, Z; 4 — operating load 
(pressure force of the liquid under discharge and suction); 5 — gravity; 6 — type of contact — connected (no gap) surfaces with 
shared nodes 


Research Results. Diagrams of equivalent von Mises stresses are shown in Fig. 6. 
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yon Mises (N/mm*2 (MPa)) 










396,392 
Max 






396.392 


363.361 






330.330 
297.299 
264.268 


231.237 





198.206 
165.175 
132.144 
99.113 
66.082 
33.051 
0.020 


— Yield strength: 325.000 





-— ae, 


c 
0.820 





c—_. _ 











b) 


Fig. 6. Diagrams obtained as a result of the strength analysis of the TWS 600 plunger pump body: a) diagram of equivalent stresses, 
b) diagram of the yield factor 


Yield factor is determined from the formula: 
O 
C=— 


O35xB 





where o, — ultimate strength of 09G2S steel, o, = 470 Mpa’; O3x3 —- Maximum equivalent stress, 


Oyen = 396.392 MPa. 


z 

5 

TN 

S 

o) 

a K — 1,18 

~ =sacsap = le 

= "396,392 

— 

© As can be seen, the equivalent stresses’ have the greatest intensity in the lower region of the pump body 
Py design. 

a= The minimum yield factor of the material is 0.82 (see Fig. 6 5). 

78 ' Zubchenko AS, Koloskov MM, Kashirsky YV, et al. Steel and Alloy Grade Guide. Moscow, 2003. 784 p. (In Russ.) 


* Feshchenko VN. Constructor reference. Book 2. Design of machines and their parts. Moscow, 2017. 400 p. (In Russ.) 
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The strength calculation of the plunger pump body was carried out for the third position of the mechanism 
(Fig. 2 5). For other positions of the mechanism (see Fig. 2 a, b) the strength calculation is performed similarly. 
The results of the strength calculation of the body for the three positions of the mechanism are shown in Table 3. 
Table 3 
Results of the strength calculation of the plunger pump body 


Mechanism 7 Equivalent stresses, Safety factor for 
- Plunger position 


1 


2 


Discussion and Conclusions. The reactions in the crankshaft supports with account for the forces generated by 





II Ls i. 284.783 1.141 1.65 


the plunger depending on its operating mode and the position of the crank are described. The forces acting on each of 
the plungers and the resulting reactions in each of the supports are determined. 

For the strength calculation of the plunger pump body made of 09G2S steel, the FEM in CAD SolidWorks 
Simulation is used. When assessing the body strength, it has been observed that in the third position of the mechanism, 
the lower region of the structure has minimum coefficients of yield factor and strength factor (see Table 3). This area 
requires strengthening since the coefficient values are insufficient for the specified operating conditions. Also, there are 
unloaded areas in the body, where the values of safety and yield factors are several times higher than the recommended 
ones. In the future, it 1s required to optimize the plunger pump body: to strengthen the areas of the structure with 
unsatisfactory coefficients of safety and yield factors, and in places where the coefficients exceed the recommended 
values, to use a metal of a smaller thickness. 

Thus, the strength calculation results can be used to optimize the design of the pump body under standard 
operating conditions. 
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Power analysis of chain transmission with gear chain and involute sprockets 


S. B. Berezhnoy, G. V. Kurapov 
Kuban State Technological University (Krasnodar, Russian Federation) 





Introduction. The power analysis of a chain transmission with a toothed chain and involute sprockets considers the 
centrifugal forces and the friction coefficients between the plate and the sprocket tooth. The work objectives are to 
determine all meshing forces, and to calculate the coupling coefficient of the gear chain with the involute sprocket in 
the drive gears. 

Materials and Methods. When evaluating the traction capacity of a chain transmission, such parameters as the power 
analysis and the coupling coefficient of the gear chain with the sprocket are important (it shows what fold the pressure 
on a given tooth is greater than the pressure perceived by the tooth in front). In the presented paper, the following 
diagrams are visualized: the arrangement of the plates in gearing with the involute sprocket teeth and the meshing 
forces. The factors that affect the involute profile of the sprocket tooth are considered. This includes the weight of the 
chain plate package and the force: centrifugal, friction, normal pressure and tension. At the same time, changes in the 
coupling coefficient for the subsequent teeth involved in traction are taken into account. The balance of the links i and i 
+ 1 of the gear chain is studied in the coordinate system XOY with the center on the axis of rotation of the involute 
sprocket. The method enables to determine all the desired forces through the geometric calculation of the values of the 
angular transmission parameters. Using the equations obtained, the following parameters are specified: the coupling 
coefficient B;, the tension of the driving branch S; and the slack branch Sy. 

Results. A patented transmission stand with a gear chain and involute sprockets is presented. The tests carried out on it 
validated the study results of a chain transmission with a toothed gear and involute sprockets with the specified 
parameters. The correctness of the power analysis of the transmission with account for the centrifugal forces and the 
friction coefficients of the plates and the sprocket teeth was proved. 

Discussion and Conclusions. It is noted that the centrifugal forces and the friction coefficients during engagement affect 
significantly the traction capacity of a transmission with a toothed chain and involute sprockets. The data obtained can 
be used to accurately estimate the traction capacity of such gears. 


Keywords: chain transmission, toothed chain, involute sprocket, chain plate, traction capacity, coupling coefficients, 
centrifugal forces, friction forces, friction coefficients, joint angular velocity, driving side, slack side. 
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Introduction. Chain transmissions with toothed chain are not sufficiently studied. Their release is limited, but 
they are increasingly used in production. Specifically, such gears are used to complete drive mechanisms in machines, 
precision machine tools, and other equipment. This paper considers the power analysis and calculation of the traction 
capacity of gears with a toothed chain and involute sprockets. 

The traction capacity determines the operability of the involute gear chain transmission, which depends on a 


number of factors’. Among them there are: 


' Kurapov GV. Performance study of a chain transmission with a gear chain and involute sprockets: Cand.Sci. (Eng.) diss. Krasnodar, 2019. 173 p. (In 
Russ.) 
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— transmission geometry, 

— limiting contact stress of the weakest element of the pair “chain plate — sprocket tooth’, 

— coupling coefficient of the toothed chain with the sprocket teeth, 

— slipping of the toothed chain along the sprocket teeth, 

— permissible specific pressure inside the joint, 

— impact resistance of gear chain elements’ [1-3]. 

Materials and Methods. The coupling coefficient is one of the main factors of the traction capacity of a chain 
transmission equipped with a toothed chain. During the operation of such a transmission, engagement occurs and forces 
act that move the chain plates along the working and occipital profiles of the sprocket teeth [4-7]. On the arc of the 
involute sprocket with a toothed chain, any of its links (link 7, i+ 1, i+ 2, etc.) experience: 

— friction force N,f of the chain on the sprocket tooth, 

— tensile strength QO; of adjacent links, 

— normal pressure force N; of the tooth sprocket profile (Fig. 1). 





Fig. 1. Arrangement of plates in engagement with the involute sprocket teeth 
In the literature, the coupling coefficient is calculated from the ratio of the tensile forces QO; and Q;+; in the 
adjacent chain links: 
pk. 
QO... 


This coefficient is not constant, so the plates of the gear chain move along the profiles of the sprocket teeth 


(1) 


throughout the entire arc of contact. Fig. 2 shows the engagement of the first package of plates (link) of the chain with 
the sprocket tooth. 


' Semenov VS. Availability assurance of roller chain drives during their operation. In: Proc. All-Russian Sci.-Pract. Conf. of postgraduates, doctoral 
students, and young scientists. Maykop, 2017. P. 121-125. (In Russ.) 


Machine building and machine science 


CO 
Oo 


http://vestnik-donstu.ru 


Advanced Engineering Research 2021. V. 21, no. 1. P. 82-88. ISSN 2687-1653 








Fig. 2. Diagram of the forces acting in the engagement of the gear chain plate with the sprocket 


The involute profile of the sprocket tooth is affected by the following factors [8-14]. 
1. Friction force N;f, which prevents the movement of the plate pack to the top of the sprocket tooth. Here, f— 
the friction coefficient. 
2. Weight of the chain plate package G, = m;g. Here, g — the acceleration of gravity, m; — chain link weight 
depending on the pitch’. 
3. Normal pressure force N; of the sprocket tooth profile concentrated at the point of engagement K;. 


4. Tensile forces Q; and Q;,; of the adjacent links of the gear chain. 


Mo 


- Here, m — angular speed of rotation 





5. Centrifugal force with account for the mass of the chain link F, = 


of the sprocket. R — radius of the central part of the link (represented as a circle that changes position when the 
sprocket rotates) on the involute profile of the sprocket tooth. 
However, when the chain moves from the driving branch to the driven branch, the tensile forces QO; and QO. 
change by the value of each value of the coupling coefficient B;; so, it needs to be determined more precisely. 
Let us build the XOY coordinate system. It passes through the center of the axis of rotation of the involute 
sprocket. The vector of the X-axis coincides with the vector of the linear velocity of the gear chain in the transmission. 
Let us analyze the equilibrium of the first link i and the second i+ 1. 
Balance of forces on the axis_X ().X = 0) will be: 
-O,+ Nyy — Nig f —Foy + ¥osay + Quay = 9. (2) 
Balance of forces on the axis Y ()}Y = 0) will be: 
G+ Go —Foyy Ny “Nyy f tO EF. 


3i 3itl a y itl Y 


=0. (3) 


At }.X = 0 and )'Y = 0, the gear chain plates will be stationary. If ).X = 0 and }:Y¢ 0, then the plates will move 
along the tooth profile. 


* GOST 13552-81. USSR State Standard. Driving toothed chains. Specifications. USSR State Standards Committee. Moscow: Standartinform; 
1987. 12 p. (in Russ.) 
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Let us determine the forces acting on the chain links. Fig. 2 shows: P — the center of mass of the gear chain 


link; d,,— the diameter of the circle of the contact point K of the package of plates with the profile of the sprocket tooth, 





d . ee t 
d. =——— ; d,— the diameter of the dividing circle of the sprocket, d, = ———. 
cos® _ 180 
sin — 
vA 
Along X-axis: Along Y-axis: 
N,, =N, -sinB, N,, =N, -cosB, 
Nf =N,f-cosB, Nf =N, f-sinB, 
Rag =i “sIn®, Es =F, "COS® , 
Eas = La ‘sing, ee a ee ‘COS®, 
Ox =0.,,, “COSO, Oey =0., ‘sIN@. 


TU 
Here, B — the angle of the tooth profile of the involute sprocket; ) = s — half of the tooth pitch angle of the involute 


sprocket. 


Values of angular parameters: 


K=OQ+Y, (4) 
u=n—y—-O , (5) 
v=~+y+0., (6) 

2 

TU 
=—-@O+0., 7 
p ‘ o+90. (7) 

TU 

—~—_g-y. 8 

B ; O-Y (8) 


Here, y — half of the tooth wedge angle of the involute sprocket. 
Y = O.,. 7 0, ° (9) 


a= arccos (10) 


0, — the angle between the axis OO;, of the tooth and the line connecting point K, to the center O of the sprocket (see 
Fig. 2). 


S 
oe 11 
7 (11) 


S;, — the thickness of the involute sprocket tooth at the point of contact K and the gear chain (see Fig. 1). d,; = mzcosa 
— the diameter of the main circle of the involute sprocket, where m — the sprocket module for the toothed chain 
transmission, z — the number of sprocket teeth. 


Substituting all the values included in the equations (2) and (3), we obtain: 


2G,—F', -cosp—N, f -sinB—N, -cosB+Q,, -sing—F’,,, -cosp=0, (12) 
(Ff, -cosp+N, -cosB+N, f-sinB+F' .., -cosp—2G, ) 
i (13) 
sin @ 
QO =N, -cosB—N, f -cosB—F’, -sing+F’,,, -sing+Q,, -cosp=0. (14) 
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Transforming the equations (13) and (14), we obtain the coupling coefficient B; of the toothed chain with a 


given involute sprocket tooth: 


O. _ N, -cosB—-N, f-cosfP-F, -sng+F ,, -sing+Q,, -cosg 
B.=— = si 9 AY  ——_). (15) 
O., Fi, -cosg+N, -cosB+N, f-sinf+F’ ,, -cosg—2G, 


If the tension value of the driving branch of the chain transmission S, is known, then the tension of the driven 


branch S,is determined with account for the overall coefficient of the chain — sprocket tooth coupling: 


S 
oe 16 
a (16) 
S 
S,=—, 17 
a} (17) 
B=B,B,-B,...B,. (18) 


Here, B, — the overall coupling coefficient of the chain and the sprocket; n — the number of teeth of the driving 
involute sprocket. 

Research Results. To verify the calculations, a stand’ protected by a patent of the Russian Federation was 
designed and manufactured (Fig. 3). Bench tests of the chain transmission under study were carried out. The 
parameters: 

1) electric motor power N,, = 12 kW, 

2) toothed chain pitch ¢,= 19.05 mm, 

3) the number of teeth of the driving and driven sprocket z, = z)= 23, 


4) speed of the driving and driven sprocket n; = n, = 1000 rpm. 





Fig. 3. Transmission test bench with toothed chain and involute sprockets: 
1— driving sprocket, 2 — driven sprocket, 3 — toothed chain, 4 — electric motor, 5 — generator 


For the bench tests, calculations were performed with account for the reference data” * * °. The results are 


presented in Table 1. 


' Berezhnoy SS, et al. Test bench for chain gears. RF Patent no. 147242, 2014. (In Russ.) 

*Chernavsky SA, et al. Design of mechanical gears: study guide for technical colleges. Moscow: Al'yans; 2008. 590 p. (In Russ.). 

* Uchaev PN, et al. Course design of machine parts based on graphic systems: study guide. Krasnodar: Izd-vo TNT; 2013. 428 p. (In Russ.) 
* Gotovtsev AA, Kotenok IP. Design of chain gears. Reference guide. Moscow: Mashinostroenie; 1982. 336 p. (In Russ.) 

> Feshchenko VN. Handbook for designers. Book 1. Machines and mechanisms. Moscow: Infra-Inzheneriya; 2016. 400 p. (In Russ.) 


Berezhnoy S. B., et. al. Power analysis of chain transmission with gear chain and involute sprockets 





Table 1 
The obtained data for determining the traction capacity of a toothed chain transmission 


Notation and dimension 
Tension of the working transmission branch 4662 
Tension of the idle transmission branch 270 


Centrifugal force of the chain jaere sl 29.82 





Total coupling coefficient of the chain and the leading B 17.2 
involute sprocket ) 


The proposed approach to determining the coupling coefficient has confirmed the correctness of the force 
calculation with account for the centrifugal forces and the friction coefficients of the plate package and the sprocket 
tooth. 

Discussion and Conclusions. In this paper, all the forces acting in the engagement of the chain plate package 
and the sprocket tooth are determined. The procedure of power calculation of a toothed chain transmission is presented. 
The dependences for determining the impact of centrifugal forces and friction coefficients on the final calculation of the 
coupling coefficient are obtained. To confirm the theory, the chain transmission was studied on a test bench. 

So, the centrifugal forces and the meshing friction coefficients significantly affect the traction capacity of the 


transmission with a toothed chain and involute sprockets. 
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On detection of crack-like welding defects by existing quality control methods 
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Introduction. The research materials devoted to the operability assurance of welded structures of hazardous technical 
devices — components of oil and gas production equipment are summarized. The factors determining the operational 
strength of welded joints and structures are systematized and diagrammed. The causes for the decrease in their service 
properties are described. We have defined the role of volumetric and crack-like welding defects, which were not 
detected during quality control at the manufacturing stage and in determining the residual life while in operation, in 
ensuring the life cycle of the structure. Measures to improve the practice of detecting crack-like defects under visual- 
and-dimensional and ultrasonic methods of quality control of welded joints are proposed. 

Materials and Methods. In accordance with SNiP 2.05.06—85, the structures of main oil and gas pipelines, vertical 
cylindrical tanks, and other oil and gas production equipment (OGPE) are made of dead-killed and semikilled low- 
carbon and low-alloy steels with tensile strength of up to 686.5 MPa. Regardless of the class and strength level of steel, 
it should be well welded by all methods prescribed by standard process documentation (SPD). At the same time, 
regardless of the state of supply, the carbon equivalent C.,, should not exceed 0.46 %. The research methods are 
calculation-experimental ones. To calculate the stress-strain state of welded joints, the following methods were used: 

— finite elements (FEM), 

— fracture mechanics using the stress intensity coefficient K,. 

Methods of mathematical statistics were used to estimate the geometric dimensions of crack-like welding defects. 
Welding defects were detected by standard quality control methods prescribed by the SPD and GOST standards. 
Results. Welded structures operating in the fatigue mode are considered. Some factors characteristic of welded joints 
are analyzed. It is shown how they affect the formation of strength performance properties. The role of dangerous 
crack-like defects, which with high probability can be formed in welded joints in the manufacture of structures and 
during operation, is established. Often, the reason is a decrease in mechanical properties due to aging and loss of 
plasticity caused by accidental mechanical actions. Note that it is impossible to identify these defects by existing quality 
control methods, both during the control process under production and during diagnostics while in operation. This 
reduces the accuracy of predicting the operational life of the welded structure. 


Discussion and Conclusion. It is proposed to include the following requirements in the SPD: 

— to the quality of welded joints of hazardous technical devices of OGPE, 

— to detection of sharp crack-like defects with a radius of curvature in the range of 0.1-0.25 mm. 

However, the reliability of detection of such defects by the control methods used remains low due to the human factor. 
To increase the reliability of detection of sharp crack-like defects, the directions of improvement of ultrasonic quality 
control of welded joints are determined. 


Keywords: welding, welds, structural strength, welding defects, non-destructive testing, visual and dimensional control, 
ultrasonic testing, validity. 
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Introduction. Scientific studies [1-7], as well as standards*, regulatory’ and specification documents’ of 
relevant organizations are devoted to the problems of determining the strength of the operated welded metal structures. 
Factors that determine the strength characteristics of weld joints are numerous and diverse. The destruction of welded 


metal structures under operation depends on: 


— operating conditions (level of calculated loads, loading scheme, joint effect of various forces, action of 


unaccounted loads, etc.); 


— structural and technological causes. 







Weld metal inhomogeneity 


Reduce the weld mechanical properties 






Weld joint structural design 
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Fig. 1. Factors determining the operational strength of weld joints and structures 


The generalization of numerous studies suggests that the load-bearing capacity of welded structures is affected 
by the following conditions. 

I. Heterogeneity of the metal structure and mechanical properties of the weld and the near-weld zone. 

II. Welding defects, mainly, crack-like ones. These are undercuts, cracks, rolls, sharp inclusions, faulty 


penetration, incomplete fusion, etc. (Fig. 2). 


* Gas-shielded arc welding. Welded joints. Main types, design elements and dimensions. Federal Agency for Technical Regulation and Metrology. 
Moscow: Standartinform; 2001. 37 p. (In Russ.) 

* RD-25.160.10-KTH-001-12. Welding technology instruction for the construction and repair of steel vertical tanks of “Transneft” JSC. 112 p. 
URL: https://www.studmed.ru/rd-25-160-10-ktn-001-12-instrukctya-po-tehnologii-svarki-pri-stroitelstve-i-remonte-stalnyh-vertikalnyh-rezervuarov- 
chast-1_5cb8180807a.html (accessed: 10.10.2020). (In Russ.) 

° Gazprom Industry Standard 2-2.4-083-2006. Regulatory documents for the design, construction and operation of Gazprom facilities. Instructions on 


non-destructive methods of quality control of weld joints during the construction and repair of field and main gas pipelines. URL: 
https://files.stroyinf.ru/Index2/1/4293831/4293831880.htm (accessed: 10.10.2020). (In Russ.) 
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Fig. 2. Defects in the butt weld joint 
III. Imperfection of the geometry of the weld cross-sections. The decrease in strength is mostly due to the 
sharpness of the interface between the base and the deposited metal (Fig. 2, 3): the smaller the transition radius r and the 
greater the angle of approach 9 from the base to the deposited metal, the lower the fatigue strength of the weld joint. 


This pattern is valid even in the absence of welding defects. 
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Fig. 3. Stress distribution in the interface zone of the base and deposited metal, at the apex of sharp defects 


The weld joint design is determined by: 

— the shape and dimensions of the butt edges of the parts; 

— relative position, configuration, length of welds; 

— type of weld joint according to GOST. 

Fig. 4 shows the stress distribution in lap joints with flank and composite welds. 
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Fig. 4. Stress distribution in lap joints: a) with flank welds; b) with composite weld 
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Table 1 shows the values of the theoretical stress concentration coefficient K, depending on the length and 
width of the tie plate: 


_ 6max 


D D 
Ku = 2 = 66-2 -cth (4,62). (1) 
Table 1 





K,, values depending on the length and width of the tie plate 





IV. Residual welding stresses (RWS). 

V. Absolute dimensions of the parts to be welded (scale factor). 

Welded joints in metal structures, even without welding defects, usually contributes to a decrease in 
operational properties. Butt weld joints are the least dangerous in this sense. Lap joints and T-joints create a stress 
concentration caused by sharp changes in the working cross-sections [8, 9]. However, for the manufacture of many 
structures, butt joints are not enough. 

Residual welding stresses can change the overall power flow distribution pattern. RWS, interacting and 
summing up with stresses from external loads, affect the performance of weld joints. This effect can also be positive, 
for example, if the vectors of RWS and stresses from external loads are multidirectional (tension + compression, etc.). 

Welded joints have the greatest impact on the fatigue of structures. Disturbing factors (Fig. 1) cause 
destruction at stresses that can be transmitted by elements of metal structures without weld joints’. 

In numerous papers devoted to the impact of welding defects on the performance of joints, recommendations 
are given to reduce harmful effects. At the same time, it is noted that some defects should be classified as harmless. For 
example, single small pores practically do not affect the fatigue strength of weld joints. In this case, the regulations on 
the presence of pores are designed to provide not the strength, but the acceptable presentation of the product. 

Traditionally, welding defects are understood as any imperfections in the geometry and discontinuities of the 
material of the welds and the heat-affected zone (HAZ). Thus, the geometry of the local interface zone of the base and 
deposited metals can also be considered a defect. 

Materials and Methods. The study of the documents regulating the presence and size of welding defects has 
shown that they practically do not require the geometry of the sections of weld joints at the junction of the weld and the 
base metal. However, it is here that destruction occurs due to the formation and growth of crack-like defects [2, 4, 10]. 
When developing standards (for example, for main pipelines operating under static and re-static loads), it is important to 
consider not only the size of the defect that weakens the working section of the weld, but also its sharpness, which 
determines the local stress concentration at the defect mouth. 

In [8], methods for predicting the stability of weld joints during fatigue work with account for the actual 
geometry of the welds and the probability of defects are presented. This approach is appropriate for the development of 
industry standards for the admissibility of welding defects considering the loading and responsibility of metal 
structures. These methods provide the determination of statistically valid characteristics of the weld joint quality. 

In metal structures under static loads, defects in welds appear if their dimensions are so significant that they 
actually weaken the working cross-sections of the joints. To determine the strength of the joints and the acceptable 
dimensions of the defects, a static strength calculation for permissible stresses is performed. 

It is advisable to approach the calculations of weld joints operating in the re-static load mode with crack-like 


defects, with regard to fracture mechanics. In this case, the stress intensity factor K, [3, 10—12] is used: 


K,=ovnl, (2) 


'Harchenko VYa, Assaulenko SS, Peredelsky VA. On the requirements for the quality of weld joints and their classification in machine-building 
structures. In: Proc. VII Sci.-Pract. Conf. on Innovative technologies in mechanical engineering and metallurgy. Rostov-on-Don: DSTU Publ. Center; 
2015. P. 341-346. (In Russ.) 
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where o — nominal voltages; / — half-length of the equivalent crack describing a typical size of crack-like defects 
defect (e.g., the undercut depth, the depth of faulty penetration, etc.). 

In this case, defects with a radius at the vertex R of less than 0.25 mm can be conventionally classified as 
crack-like defects. These are the defects for which R > 0.25 mm can be conventionally considered blunt. 

Research Results. One of the factors that determine the operational properties of welded structures are crack- 
like defects, such as: faulty penetrations, undercuts, inter-roll incomplete fusions, microcracks. 

During visual and measuring control, it is very difficult to detect sharp defects with small measured values 
(0.1-0.25 mm). This is due to the low sensitivity of the measuring instrument and the human factor since the 
measurements are performed using the controller in manual mode [13]. 

Radiography and ultrasound are the two main non-destructive physical testing methods for detecting internal 
defects. Under studying the particularly critical structures, these methods often duplicate each other since they have 
different sensitivity and reliability. To bring out sharp, undetectable defects (e.g., cracks), an ultrasound control’ is 
preferred [14, 15]. However, it is often done manually. Accordingly, the results may be affected by the human factor, so 
this method is not sufficiently reliable. 

Major challenge in assessing the product quality may arise for specialists already at the preparation stage, when 
choosing the method of ultrasonic testing (e.g., when using the mirror-shadow method, the echo technique, or the 
shadow method). To choose, you need to know the methods of operation, as well as to consider the characteristics of the 
converters in accordance with GOST R 55725-2013*. Such knowledge is not always demonstrated even by third-level 
specialists, which is revealed during certification. 

Equipment and quality control specialists should be able to correctly configure the flaw detector according to 
standard samples (SS-1, SS-2, SS-3) and according to the enterprise standard samples (ESS). Such templates are made 
of a material similar to the material of the object under control. The samples should correspond to the object in 
thickness and (if required) in radius (when checking small diameter pipes). However, the practice of conducting 
certification activities indicates that in some cases, specialists are not able to perform this work efficiently. 

During the control process, the speed and step of scanning the object are also important. If the scanning speed 
is too high, the operator risks missing the defect or inaccurately determining its coordinates, the depth of occurrence and 
the boundaries of the discontinuity. The scanning speed standards are specified in the regulatory documentation. 
However, in the field, the control procedure is complicated by temperature difference, which depend on the time of 
year, venue, etc. For example, depending on the environment temperature, the temperature of the object under control 
changes, and this can affect the viscosity of the contact liquid. On the other hand, at low temperature, heating of the 
controlled area is required. This changes some of the physical properties of the object under control and, accordingly, 
the passage of ultrasound in the product, etc. 

The main problem in summing up the results of the control and issuing conclusions is a free and subjective 
interpretation of the data obtained on the number of defects, depth and area of their occurrence. Much depends on the 
experience of the operator, his ability to distinguish the signals of defects from false or superficial ones. 

Discussion and Conclusions. It is required to restrict the influence of the human factor to increase the 
reliability of the results of visual-measuring and ultrasonic quality control of weld joints. To achieve this goal, we 
consider it appropriate to organize additional training in the process of periodic mandatory assessment of the 
qualifications and certification of control specialists. Teachers should be experts of welding production. In short-term 
training programs, it is required to summarize information on the impact of sharp crack-like welding defects on the 


strength characteristics of welded structures and their operational life. 


' GOST 55724-2013. Non-destructive testing. Welded joints. Ultrasonic methods. Federal Agency for Technical Regulation and Metrology. 
Moscow: Standartinform; 2019. 28 p. (In Russ.) 

* GOST P 55725-2013. Non-destructive testing. Piezoelectric ultrasonic transducers. General technical requirements. Federal Agency for 
Technical Regulation and Metrology. Moscow: Standartinform; 2019. 16 p. (In Russ.) 
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During certification, it is required to use special simulators for ultrasonic testing and samples of weld joints of 
structural elements with artificially created defects that simulate their real size, shape and location in test samples. This 


will provide obtaining an objective assessment of the reliability and informativeness of ultrasound control. 
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Introduction. In all types of digital communication, error control coding techniques are used. Many digital 
communication standards, such as W1-Fi and 5G, use low density parity check (LDPC) codes. These codes are popular 
because they provide building encoders and decoders with low computational complexity. This work objective is to 
increase the error correcting capability of the well-known bit-flipping decoder (BF) of LDPC-codes. For this purpose, a 
modification of the decoder is built, which enables to dynamically control one of its main parameters whose choice 
affects significantly the quality of decoding. 

Materials and Methods. The well-known bit-flipping decoder of binary LDPC-codes is considered. This decoder has 
several parameters that are not rigidly bound with the code parameters. The dependence of the decoding quality on the 
selection of the output parameters of the bit-flipping decoder was investigated through simulation modeling. It is shown 
that the decoding results in this case are significantly affected by the input parameter of the decoder — threshold 7. A 
modification of the BF-decoder of binary LDPC-codes has been developed, in which it is proposed to set the threshold 
dynamically during the execution of the algorithm depending on the error rate. A comparative analysis of the error- 
correcting capability of decoders 1s carried out by the simulation modeling method. 

Results. A lemma on the maximum value of the decoder threshold T is formulated and proved. Upper bounds for the 
number of operations are found for the original and modified decoders. A simulation model that implements a digital 
noise-immune communication channel has been built. In the model, the initial data is encoded with a given LDPC-code, 
then it is made noisy by additive uniformly distributed errors, and thereafter, it is decoded in turn by the bit-flipping 
algorithm with different threshold T parameters, as well as by a modified decoder. Based on the input and output data, 
the correction capacity of the decoders used is estimated. Experiments have shown that the error-correcting capability of 
the modified decoder in the range of the real error rate is higher than that of the original decoder, regardless of the 
selection of its parameters. 

Discussion and Conclusions. The lemma, proved in the paper, sets the upper bound on the threshold value in the 
original decoder, which simplifies its adjustment. The developed modification of the decoder has a better error- 
correcting capability compared to the original decoder. Nevertheless, the complexity of the modification is slightly 
increased compared to the original algorithm. It has been pointed out that the decoding quality of a modified decoder 


develops with a decrease in the number of cycles in the Tanner graph and an increase in the length of the code. 


Keywords: LDPC-codes, error-correcting capability, dynamic threshold, binary symmetric channel, experimental 


research. 
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Introduction. In 1963, in [1], R. Gallager first described a class of linear block codes whose check matrix 
contains a small number of nonzero elements. Such codes are commonly referred to as low-density parity check codes, 
or LDPC codes. For them, it is possible to build encoders and decoders with low computational complexity. Thus, when 
using LDPC codes, the data transfer rate is not significantly limited. Many modern studies are devoted to LDPC codes 
and their decoders [2-5]. LDPC codes are widely used in various digital communication standards, such as W1-F1, 5G, 
and optical communication [6, 7]. However, despite the popularity of these codes, some of the problems associated with 
them require research and solution. One of them is building new decoders and improving the existing ones. 

This work objective is to increase the error-correcting capability of the well-known bit-flipping decoder of 
LDPC codes (hereinafter referred to as the BF decoder). To do this, a modification of the decoder is built, which 
enables to dynamically control one of its key parameters, whose selection affects significantly the quality of decoding. 

Materials and Methods. The key parameters of binary LDPC codes are length N, dimension K and minimum 
code distance d. The information words [N, K, d] of the C-code are vectors m = (m,,mMj,..., Mx) € Fs , Where F, is the 
Galois field of cardinality 2, and the codewords are vectors € = (C1,Cp,...,Cy) € Fs’ [8]. It is convenient to set the 
LDPC codes with the check (N — K) X N matrix H. Most of its elements are zero [1], so it is more convenient to store 
it not entirely, but storing only the positions of nonzero elements rowwise. 

There are regular [9] and irregular [10] LDPC codes. In regular codes, all rows and columns of the check 
matrices contain a fixed number of single elements (k and j, respectively), otherwise the code is called irregular. For 
convenience, check matrices of regular LDPC codes will be called regular matrices, and irregular LDPC codes — 
irregular. 

Regular LDPC codes have a number of advantages: easily evaluated code parameters, easy storage of matrices, 
low computational complexity of encoding and decoding algorithms, etc. In addition, regular code decoders correct 
errors evenly, unlike irregular ones, which correct errors in some parts of the codeword worse than in others. However, 
the problem of generating regular matrices with given properties is complex, and brute-force methods are often used to 
solve it. 

To discuss the properties of the matrix H it is convenient to use the corresponding Tanner graph G = (V,E), 
where E — a set of edges, and V = SUR — a set of vertices, S — a set of rows of the matrix H, and R — a set of its 
columns [11]. Each nonzero element H, standing in the i-th row and the j-th column, defines an edge connecting the i- 
th vertex of the set S and the j-th vertex of the set R. Fig. 1 shows an example of a regular check matrix 3x6 with 


parameters k = 4 and j = 2, and the corresponding Tanner graph. 





Fig. 1. The cycle in Tanner graph and in the check matrix 


The top row of the graph vertices corresponds to the columns of the matrix H, and the bottom row is connected 
to the rows of H. An important characteristic of the check matrix H of the LDPC code 1s the presence and type of cycles 
in the corresponding Tanner graph. A cycle is a sequence of adjacent vertices of a graph in which the first and last 
vertices coincide. The length of this sequence is called the cycle length. The minimum cycle length in a graph is called 
the girth. If the graph contains no cycles, its girth is assumed to be infinite. An example of a cycle of length 4 is 
highlighted in bold lines in the graph (Fig. 1). 
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The error-correcting capabilities depend not only on the key parameters of the LDPC codes, but also on the 
structure of the check matrix H. On the one hand, the presence of cycles of small lengths (such as 4 and 6) impairs 
noticeably the error-correcting capability of the decoder [12]. On the other hand, the code that corresponds to Tanner 
graph without cycles does not correct errors, since its minimum code distance is 2. Thus, the task of constructing check 
matrices of regular LDPC codes is multiparametric. When solving it, you need to monitor the key parameters of the 
code, as well as the cycles in Tanner graph corresponding to the check matrix. 

Consider the well-known BF-decoder of the LDPC code C in a convenient form [13]. 

Input: LDPC code C with parameters [N, K,d], given by the check matrix 


hi hi ann hin 
Aw-k)1 Aw-Kj2 = Aw-Kw 
Vector ¢ =€+é@, CE C(c F’), @(€ Fs’) — error vector; p — the number of iterations of the algorithm; 7 — 


threshold value. 

Output: code vector € € C(c F¥). 

Step 1. Let the counter r be equal to zero. 

Step 2. Calculate the syndrome 5 = ¢’H’. If 5 = 0 orr = p, then go to step S. 

Step 3. Select the unit coordinates from the vector 5 = (s1, 59, ...,) Sy_x), L¢., S; = 1,i = 1,(N — K). Compose 
the set L = {ils; = 1}. Calculate h’ = (hj, h}, ..., hy), where 

hy = dier hu. (2) 

The values h,,, 1 = 1,..., N should be assumed to be nonnegative integers. Thus, h’ € NW’, where Ny = N U {0}. 

Step 4. In the vector h’ = (hj,hj,..., hy), we find all the elements h; > T. Among them, we select random h; 
and invert the bit c; of the vector c. Add a unit to the counter r and go to step 2. 

Step 5.c := Cc’. 

The research carried out in this work allows us to make some observations on the BF-decoder. 

Observation |. The input parameter p sets the maximum number of iterations of the algorithm from the 2nd to 
the 4th steps, but the decoder can recover the codeword in fewer iterations. 

Observation 2. When selecting the parameter T, the following considerations should be taken into account. If 
the parameter d of the used [N, K,d]|-code C is known, then it can be applied to calculate t — the number of reliably 


recoverable errors, and then the number of decoder iterations is limited to this value: 


p=t=|—. (3) 


2 
Here, |x] — rounding the number x to a smaller integer. If the parameter d is unknown, then it can be estimated using 


the Singleton bound [5] 
d<N-K+1 
and, using (3), we obtain 
N-K 
r= ES] 
Observation 3. The structure of the decoder is such that the recovery of the correct codeword is not guaranteed, 
even if in the noisy word c’ = € + €, no more than t errors occurred (3). 


Observation 4. In the literature, for regular check matrices in the BF decoder, it 1s recommended to select the 
threshold T depending on the weight j of the column of the matrix H, namely, T = - For irregular matrices, such 


recommendations are not given in the literature. The error-correcting capability of the BF decoder can be worsened by 
an unsuccessful selection of threshold T. If its value is large, at step 4 of the decoder, the vector h’ may not have a 
coordinate that exceeds the threshold T, therefore, the erroneous bits will not be corrected. If you select a small value of 
T in step 4 of the BF decoder, several coordinates, whose value exceeds the threshold, may appear in the vector h’. 
Among them, there may be coordinates that do not contain an error. Thus, the selection of the parameter 7 can 


significantly affect the decoding quality. 
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Research Results. We formulate and prove a lemma on the maximum possible value of the threshold T. Then 
we modify the BF decoder so that the threshold is set dynamically during the decoding process, and conduct a 
comparative analysis of the original and modified decoding algorithms. 

Lemma. Let the binary [N,K,d]|-code C be given by the check matrix H having a fixed number of 7 unit 
elements in each column. Then the maximum threshold value T for the BF decoder of such LDPC code C cannot be 
greater than 

fey a (4) 

Proof. Let the vector c’ = ¢ + € be obtained from the transmission channel, where c € C — is the correct 
codeword, @ € Fj’ — the error vector with the Hamming weight w(é). If w(é) = 0, then, in step 2, the vector- 
syndrome 5 = 0. Hence, the algorithm will go to step 5 and return é’ as the answer. In this case, the threshold value is 
not used. If w(é) > 0, then the regularity of H implies the validity of the inequality h; < j, where h;, — the elements of 
the vector h’. The inverting of the bit c, of the vector ¢’ occurs in the algorithm only if h; > T. Therefore, 

T<h, Sj. 

Thus, the formula (4) is correct. 

We will make changes to the BF decoder that will allow us to determine the threshold value dynamically, 
depending on the degree of damage to the code vector in the transmission channel. 

Input: [N, K, d]-code C given by the above check matrix (1). Vector ¢’ = € + &, where € € C(c FY), @ (€ F,) 
— the error vector; p — the number of iterations of the algorithm; T — some threshold value selected in advance. 

Output: code vector € € C(c FW’). 

Step 1. Let the counter r be equal to zero. 

Step 2. Calculate the syndrome 5 = ¢’H’. If 5 = (0,...,0) or r = p, then go to step 7. 

Step 3. Select the unit coordinates from the vector 5 = (s1, 59, ...,Sy_x), i-¢., S; = 1,i = 1,(N — K). Compose 
the set L = {ils; = 1}. Calculate h' = (hj,h},..., hy), where h; is the same as in the original decoder (2). When 
summing the value h,, | = 1,..., N, we should assume nonnegative integers. Thus, h’ € NX, rae N, = N u {0}. 

Step 4. Initialize the threshold value T = max(h;))<1..y — 1. 

Step 5. If T = 0 

Select an arbitrary element hj, 'of the vector h’ — such as hi, > T. 
Invert the bit cg. 

Step 6. Add a unit to the counter r and go to step 2. 

Step 7.¢ = Cc’. 

Observation 5. The modified algorithm generally performs fewer iterations than the BF decoder since the 
threshold is selected dynamically in step 4. Therefore, the decoder does not perform useless iterations that do not 
change the bits of the vector c’. The threshold value in the modified decoder depends on the number of errors that 
damaged the codeword, and is immediately set so that the noisy codeword c’ is guaranteed to be changed. 

Let us estimate from above the number of addition, comparison and assignment operations in both decoders. In 
the original BF decoder of the [N,K,d|-code C, p(kK +(N—K)N +1) addition operations, p(3N — 2K + 2) 
comparison operations and p((N —K)(k+1)+2N+ 3) + 1 assignment operations are performed. In the BF decoder 
with dynamic threshold, p(kK + (N — K)N +3) addition operations, p(SN — 2K + 3) comparison operations and 
p((N —K)(k+1)+2N+ 4) + 1 assignment operations are performed. Here, p — the decoder parameter that sets the 
maximum number of operations, k — weight of the code check matrix rows. Note that when implementing the 
algorithm, the multiplication and division operations are not actually used, as long as at the second step, it is convenient 
to use addition operations instead of multiplication to calculate the syndrome S. Recall that the matrix H has a sparse 
structure, and its rows are conveniently stored as a list of nonzero element numbers. Therefore, instead of multiplying 
the vector C’ by the matrix H, it is required to sum the coordinates of the vector C’, whose numbers coincide with the 


numbers of nonzero elements in the corresponding row of the matrix H. 
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Compared to the original algorithm, the modified BF decoder performs more operations, but moderately: the 
number of comparison operations has increased by p(2N + 1), assignment operations — by p, addition operations — 
by 2p. 

For a comparative study of the error-correcting capability of the original and modified decoding algorithms, a 
software tool has been created that implements a simulation model of a binary symmetric perfectly synchronized noise- 
immune communication channel according to [14-16]. To provide noise immunity, the model uses LDPC codes and BF 
decoders (original and with dynamic threshold). Errors in the channel were modeled as independent and uniformly 
distributed. 

The experiments used purposely found check matrices that specify LDPC codes. We describe the key 

parameters of these matrices using the standard notation of the key parameters of the code, as well as: j and k — the 
weight of each column and the weight of each row of the check matrix, respectively; w,,0, — 4 and 6 cycles in Tanner 
graph corresponding to the check matrix. 
Regular matrix H,;:N = 20,K =5,j =3,k =4,d =6,0, =0, 0, = 41. 
Regular matrix Hj: N = 28, K =7,j =3,k =4,d =6,w, = 0, W, = 42. 
Regular matrix H3:N = 28, K =7,j =3,k =4,d = 6,0, =0,W, = 29. 
Irregular matrix H,: N = 32, K =5,j =3,d=12,0, =0,0, = 0. 

Using these matrices, LDPC codes were constructed and simulation experiments were conducted. Fig. 2—5 
show graphs of the dependence of the error-correcting capability of the constructed LDPC codes on the error probability 
in the channel. For the rationale for selecting the threshold values T = 1 and T = 2 in the BF decoder, see Observations 


3, 4 and Lemma. 
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Fig. 2. Graph of the decoder error-correcting capability for LDPC codes given by the matrix H, 
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Fig. 3. Graph of the decoder error-correcting capability for LDPC codes given by the matrix H3 
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Fig. 4. Graph of the decoder error-correcting capability for LDPC codes given by the matrix H, 
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Fig. 5. Graph of the decoder error-correcting capability for LDPC codes given by the matrix H, 


In the range of the real error level [8, 13, 14] in Fig. 2-4, it can be observed that the BF decoder at the 
threshold value T = 2 shows better results than at T = 1, and the modified BF decoder has a better error-correcting 
capability compared to the original one. 

The decoders show similar efficiency at small values of the code length, but when it is increased, the modified 
decoder shows better results. Specifically, if the error probability in the non-noise-immune channel is 0.05, the 
difference in the error probability in the noise-immune channel when using a BF decoder with the threshold T = 2 and 
T = 1, 1s from 0.005 to 0.03 in favor of using a larger threshold value. If a BF decoder with the threshold T = 2 and a 
modified decoder are used, this difference varies depending on the LDPC code in the range from 0.001 to 0.003. If the 
error probability in the non-noise-immune channel is 0.1, the error probability in the noise-immune channel when using 
a BF decoder with the threshold T = 2 is less than with the threshold T = 1 by the value from 0.001 to 0.02. When 
using a BF decoder with the threshold T = 2 and a modified decoder, this difference varies in the range from 0.002 to 
0.01depending on the LDPC code. 

Both decoders are sensitive to the number of cycles in Tanner graph corresponding to the LDPC code check 
matrix. The greater the ratio of the number of cycles to the total number of elements in the matrix, the worse any BF 
decoder corrects errors. During the experiments, it was interesting to find out whether it is possible to increase the 
number of cycles in the matrix so that the modified decoder will show worse results compared to the BF decoder. 
Experimentally, the matrix H, containing 41 cycles of length 6 was found. The results of the study on the error- 
correcting capability of decoders for this matrix are shown in Fig. 5. Note, however, that the matrix H, contains even 
more cycles of length 6, namely, 42. The fundamental difference between the matrices H, and H, is in the density of 
units: 

— in H, — 60 unit elements per 300 matrix elements, 
— in H, — 84 units per 588 matrix elements. 

Recall that the feature of LDPC codes is the sparse structure of the check matrix, so H, 1s more typical for 
LDPC codes. 

Discussion and Conclusions. The paper considers a bit-flipping decoder for binary LDPC codes. 


Recommendations on the selection of such input parameters of the decoder as the threshold and the number of iterations 
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of the algorithm are given. A lemma on the maximum value of the decoder threshold is formulated and proved. A 
modification of the BF decoder of binary LDPC codes has been developed, in which it is proposed to set the threshold 
dynamically during the execution of the algorithm depending on the resulting syndrome. For the original and modified 
decoders, upper estimates of the number of operations are found. These estimates show that the modification 
complicates the decoder only slightly. The conducted simulation experiments demonstrate better error-correcting 
capability of the modified decoder in relation to the original one. The experiments also showed the dependence of the 
decoding quality on the degree of matrix sparsity and the number of cycles of length 6 in the Tanner graph 
corresponding to the check matrix of the LDPC code. Thus, the problem of constructing check matrices with a small 


number of short cycles arises, which is the subject of further research. 
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Introduction. This paper is devoted to the development of a method for estimating the current time-frequency 
parameters of each of a set of simultaneously and independently operating generators in the radio electronic system. A 
general case is considered, in which the deviation of true values of the nominal generator parameters from the assumed 
values of these parameters is determined not only by random, but also by long-term frequency deviation. The work 
objective is to generalize the method for estimating the time-frequency parameters of signals (long-term nominal 
frequency and current frequency deviation from the nominal value) based on the simultaneous measurement of the 
phases of signals generated in the system of independently functioning generators. The research task is to consider a 
system of simultaneously and independently operating generators. Each of them generates harmonic signals, whose 
time-frequency parameters, such as the average frequency, are constant during a certain interval of observation. But 
herewith, these time-frequency parameters are known with insufficient accuracy due to the influence of external factors 
(changes in temperature, pressure, supply voltage, etc.). It is required to obtain estimates of the time-frequency 
parameters of signals (the duration of the measuring interval, values of the long-term frequency and the standard 
deviation) from the results of measurements of the phases of signals formed by the generators at measuring intervals 
belonging to the observation interval, within which the average frequency remains constant. 

Materials and Methods. A system of simultaneously and independently functioning generators is considered. The long- 
term value of signal frequency for each of the generators over the observation interval remains constant, but it is known 
with some margin of error. During the observation interval, several measurements of the signal phase of each of the 
generators are performed. At the same time, the current values of the signal frequency and the duration of the measuring 
interval have random deviations from the long-term values, and follow the normal distribution law with zero 
mathematical expectation and a known variance. The estimation of time-frequency parameters based on the results of 
measuring the signal phases is carried out using a multidimensional likelihood function. The maximum is found on the 
base of solving the redefined system of linear algebraic equations. 

Results. A new mathematical model and a numerical-analytical method for determining the time-frequency parameters 
of signals are developed. They take into account both the long-term constant frequency deviation and short-term 
random deviations. 

Discussion and Conclusions. The results obtained can be used under the development and creation of data-measuring 
and information-telecommunication systems, including geographically distributed systems. The resulting estimates of 
the time-frequency parameters enable to increase the signal frequency stability and, accordingly, to improve the 
accuracy of measurements and the quality of information transfer. 
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Introduction. The constant growth of requirements for audio and video information transmission, for data 
formation and transmission systems, necessitates an increase in stability of the time-frequency parameters of signals 
generating in radio-electronic systems (RES) [1]. The information-telecommunication systems and data-measuring 
ones are typical examples of such systems. One of the components in providing high efficiency of such systems 1s 
associated with the formation of signals with highly stable time-frequency parameters. The latter is relevant, in 
particular, for the systems of information transmission with complex signals, radar and radio navigation systems, 
and audio and video information transmission systems as well. 

Currently, the main approach in the creation of RES is a modular construction principle, which determines the 
use of functionally completed blocks and devices connected to each other. The application of this approach leads to 
presence of a large number of high-frequency signal (HF signal) generators in any radio-electronic system. At the same 
time, functional completeness of each of the devices included in the system enables to consider the HF-signal 
generators of these devices as units operating simultaneously and independently on a certain time interval. 

Despite the fact that signal generators in various devices have different parameters, the presence of a large 
number of simultaneously and independently functioning generators in the RES provides determining the current values 
of the time-frequency parameters of these signals through measuring and subsequent processing of the generated signals 
phases’ [2-6]. The latter makes it possible to estimate the current values of the time-frequency parameters with higher 
accuracy. In turn, the obtained values of the time-frequency signal parameters are associated with the parameters of the 
generators that are forming them. That allows either to stabilize the generators frequency, or to consider its current 
value during the subsequent signal processing [2-6]. 

Well-known solutions to the problem of estimating the current time-frequency parameters take into account, as 
a rule, only random components of the frequency deviation due to the influence of various factors. In this case, the 
values of such parameter as the nominal (mean) frequency of the generator are considered known. At the same time, in 
many cases, long-term frequency deviations associated with both the influence of external factors and the technology of 
generator production are not analyzed in the papers under examination. 

Thus, the solution to the problem of estimating the time-frequency parameters of signals, such as the nominal 
(mean) frequency and random deviations from the average frequency of generators, is a challenge. 

The paper objective is the generalization of the method of estimating the time-frequency signal parameters (the 
long-term nominal frequency and the current frequency deviation from the nominal value) based on the simultaneous 
measurement of the signal phases forming in the system of independently functioning generators. 

As in the paper [3], let us consider the system of N +1 generators, each of which forms the harmonic signals, 


whose time-frequency parameters, such as average frequency o,, are constant, but known with insufficient accuracy 


due to the influence of external factors (changes in temperature, pressure, supply voltage, etc.) during a certain 
observation interval. 

It is required to obtain estimates of the time-frequency parameters of signals (the duration of the measuring 
interval, frequency values and RMS (root-mean-square, or standard, deviation) of signals) from the results of 
measurements of the phases of the signals formed by the generators at the measuring intervals belonging to the 
observation interval, within which the average frequency values remain constant. 

Mathematical model and solution method. For the system of N+1 generators under consideration, let us 


present the signals formed by each of the generators, as 
s,(t) =A, cos(@, -t+9,), n=1L..,N41. (1) 
Using one of the generators (hereinafter the N+1 generator), we will form measuring intervals of nominal 
duration ¢°’ (m=1,...,M ), such that during all the specified intervals, the average frequency and RMS of the 


frequency of each of the generators can be considered constant. 


We represent the measured values of signals phase ®,, of each of N generators (n =1,...,N ) on the m-th 


measuring interval (m=1.,...,M ) as follows: 


* Gabrielyan DD, Prygunov AG, KhutortsevVV, et al. Generator frequency stabilization method. RF Patent no. 2219654, 2003. (In Russ.) 
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®,,, =(@, +Ao,,, (1 +.At,), n=1..,N, m=1..M, (2) 


n, 


where Aw, — is random deviation of the n-th generator frequency on the m-th measuring interval from the 


average value w ; f°” and Af, are, respectively, nominal duration and deviation from the nominal value of the m-th 


measuring interval duration. 
For the mean value frequency, taking into account the modern production technology of HF—generators, we 
can write 
o =o” +60, 60 Kk”, n=1,...,.N+4+1, (3) 


(0) 


where the values @.’ are known, and 6@, , with account for the influence of external factors, are unknown. 


Taking into account the expressions (3), the relation (2) takes the form 


®,,, =(@ +80, +Aa,,, (1 +At,), n=1,...N, m=1..M. (4) 


nm 


After linearization (the terms 6@,-Az, and Aw, -At,, are discarded), the expression (4) enables to estimate 


the random deviation of the frequency of the n-th generator on the m-th measurement interval as follows: 


Aw, ,, =—2 #4, nahN, m=l,...M, (5) 
t 


es 


where ®' es 


=o. 
The frequency deviations Aw,,, (1=1,...,.N,m=1,...,M ) and deviations of measurement interval duration 


At, (m=1,...,M ) follow the normal distribution law [7]: 


4 (Ao,,,.) 1 (At,,) 
P(e) Brae | aa -oy |)” ae, | 30,0) | » 


where o” are known values of the relative instability of the frequency of the n-th generator (7 =1.,..., N ). 


The given relations describe the mathematical model of the system of generators that function simultaneously 
and independently. 

Solution method. Taking into account the relations (5) and (6), we can write the multidimensional logarithmic 
likelihood function: 


n,m 


OP —O —8@ -1 —@™-At ) 
L (dq, At) = Sy infPa+ ing, + bean Pam ~ 80, “te — 0) “Aty ) 


2 
m=l n=l ) (0 . wo” ‘ ie ) 


(7) 


which includes the vectors 0@ and At, whose elements are unknown values, respectively, 60, (n=1,...,N ) 
and At (m=1,...,M ). 
The estimates 6, and At, are found from the maximum condition (7) and correspond to the solution of the 


system of linear algebraic equations (SLAE). 


AL (de, At) 
————=(0, n=],...,N, 
00, (8) 
AL (de, At) 
———§_—=0, m=].,...,M. 
OAt, 
The equations in (8) have the form 
(D,,, —B° Sa, 1 — a - At, )- 
nym nim n m - m m = , n = 1,...,N : 
(a y ) 
(9) 
(® —@ ~§@ -¢% —@ - At )-@ (0) 
ee ee =0, m=l,...,.M 


(0. -0”-4”) 


In a matrix form, the SLAE (9) can be presented as: 


A,, A, \{ do B, 
= ; (10) 
Ag 2X55 J AL B, 
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2 ; : 
where 4. is the block with dimensions N-MxN with elements a“ =(o0°-@°) ; 4, is the block with 
1,1 n n 1,2 


nm,m 
~2 -1 

, , ’ 12 0 0) (0 
dimensions N- MxM with elements a‘? =(o0) -(@ -1) 


nmm 


T 
A,, is the block with dimensions N-M x N with elements a‘, = ioe) 


—2 


A,, is the block with dimensions N- MxM with elements a‘? = (o(” t)”) : 


nmm 


©, - 0” 
B, is the block with dimensions V-M x1 with elements 5°? =—_—"—__*"___; 
6) .) 4 
©, - 0” 


B, is the block with dimensions N-M x1 with elements b‘“’? =——"_*"_ ; 
(0 a) ; Ae 
T is the sign of the matrix transposition. 


The system of equations (9) contains N of the unknown components 6,and M of the unknown 
elements At, . The total number of measured values of signal phase is N-M . The estimation of their values, taking into 


account the measurement errors, should be carried out by the least-square method (LSM) [8-10]. In this case, the 


condition N-M >N+M must be met, and the number of measurement intervals must satisfy the next condition: the 
number of measurements M must fulfil the condition M > N/(N-1). 


The representation of SLAE using the expressions (9) and (10) defines the numerical-analytical method for the 
problem solving. All elements of the matrix have an analytical representation. At the same time, when passing to the 
system of normal equations, as a rule, used in the LSM, obtaining the analytical expressions is also not difficult. 


However, the inversion of large-dimensional matrix can be performed only using the numerical methods [11-13]. 
The estimates (do, ) and (Az, y obtained from the solution to the system of equations, represented by the 


expressions (9), allow us to determine current time-frequency parameters of the generators and signals forming by 
them. 

Research Results. We propose a mathematical model describing a simultaneous and independent functioning 
of the generator system, and a numerical-analytical method for determining the time-frequency parameters of the 
signals with account for both the long-term constant frequency deviation and the short-term deviation of random nature. 
This approach enables: 

— to evaluate, according to the results of measuring the phases of signals formed by simultaneously and 
independently functioning generators, not only random deviations in the frequency of the generated signals, and to 
obtain the estimates of the average long-term frequency of each of the generators; 

— to exclude an instability influence of the time interval duration of measurements on the resulting 
estimates of the signal time-frequency parameters [14-16]. 

Discussion and Conclusions. The results obtained can be used under the development and creation of data- 
measuring and information-telecommunication systems, including the geographically distributed systems. The resulting 
estimates of the time-frequency parameters provide increasing the signal frequency stability and, accordingly, 


improving the accuracy of the measurements and the quality of information transmission. 
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